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ABSTRACT 


An indirect adaptive control algorithm has been studied 
to control physical systems with parameter uncertainties. 
Although the particular algorithm investigated is applicable 
to a wide class of discrete time linear systems, parameter 
convergence and therefore global stability of the whole system 
is guaranteed only if the external excitation contains a 
sufficiently large number of frequency components. 

This research study has also investigated the possibility 
of stopping the identification procedure when the parameter 
error becomes sufficiently small, so the controller auto- 
matically turns itself from adaptive to time invariant, while 
still guaranteeing global stability of the closed-loop system. 
In this way the adaptive controller might be activated only 
when its gains do not provide satisfactory performances. 

Computer programs of the indirect adaptive control have 
been written for simulation purposes using both recursive 


least squares and projection algorithms. 


TABLE OF CONTENTS 


ie tele ht) ————{——{——{—[—{[——{$——$—————— a 
A. GAIN SCHEDULING ---------------------------- 
B. MODEL-REFERENCE ADAPTIVE SYSTEMS ----------- 
C. SELF-TUNING REGULATORS AND POLE PLACEMENT -- 
ia ZEROES OF SAMPLED DATA SYSTEMS ----------------- 


oe COME BOMATN ANALYSTS “OF THE ZEROES OF 
SAMPLED DATA TRANSFER FUNCTIONS ------------ 


B. FREQUENCY DOMAIN ANALYSIS OF THE ZEROES OF 
SAMPLED DATA SYSTEMS ----------------------- 


III. INDIRECT ADAPTIVE CONTROL FOR DISCRETE TIME 
SYSTEMS ---------------------------------------- 


A. POLE PLACEMENT DESIGN BASED ON INPUT- 
OUTPUT MODELS ------------------------------ 


B. DIOPHANTINE EQUATION ----------------------- 
C. PARAMETER ESTIMATION ----------------------- 
D. PERSISTENCY OF EXCITATION ------------------ 
Peo bee Mere Ro hotENCyY OF EXCTTATION ==-=-- 
iv. SIMULATION STUDIES ----------------------------- 
APPENDIX A: DESCRIPTION OF THE COMPUTER PROGRAMS ---- 
APPENDIX B: COMPUTER PROGRAM CONDIS ----------------- 
APPENDIX C: COMPUTER PROGRAM RCIOR ------------------ 
et Dee ee COMPUTER PROGRAM REMOP™ —=——>—>-—-—---=Se=- 
LIST OF REFERENCES ----------------------------------- 


INITIAL DISTRIBUTION LIST ---------------------------- 


10 


14 


14 


22 


28 


30 


26 


oo) 


44 


46 


50 


el 


83 
eel 


I. INTRODUCTION 


One of the most challenging, interesting and active fields 
of Automatic Control is Adaptive Control. To implement high- 
performance control systems when the plant dynamics are poorly 
known Or when large and unpredictable variations occur, the 
control engineers prefer to use an important class of control 
systems called Adaptive Control systems. Adaptive Control 
comes from a desire and need for improving performance of 
complex engineering systems with large uncertainties. It is 
especially important in systems with many unknown parameters 
that are changing with time. Also, it can be defined as a 
special type of nonlinear feedback control, as a nonlinear, 
nonautonomous dynamic system. 

An adaptive controller can change its behavior in response 
to changes in the dynamics of the plant and the disturbances. 

The term adaptive control has been used at least from the 
beginning of the 1950's. With recent advances in micropEees 
essor technology, it has become feasible to implement adaptive 
algorithms efficiently in real time at a reasonable cost. 

There are three schemes for parameter adaptive control: 
gain scheduling, mogel reterence control and self-tuning 
regulators. The starting point is an ordinary feedback ¥ceus 
trol loop with a process and regulator with adjustable 


parameters. But the main problem is to find a convenient way 


of changing the regulator parameters in response to changes 
in. process and disturbance dynamics. The following schemes 
differ only in the way the parameters of the regulator are 


adjusted. 


A. GAIN SCHEDULING 

Gain scheduling depends on finding auxiliary process 
variables correlated with the changes in plant dynamics. [In 
this way it is possible to reduce the effects of parameter 
variations by changing the parameters of the regulator as 
functions of the auxiliary variables. 


The block diagram of this scheme is given in Figure 1.1. 
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Figure 1.1. Block Diagram of Gain-Scheduling System 


B. MODEL-REFPERENCE ADAPTIVE we donors 

In this type of adaptive system, a reference model speci- 
fies the desired performance, and tells how the process output 
should respond to the command signal. A block diagram of 
model reference system is given in Figure 1.2. As seen, from 


this figure, the reference model is part of the control system. 
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Figure 1.2. Block Diagram of the Model Rerer-ice 
Adaptive System 


This adaptive system has two loops: the lower loop con- 
tains Regulator and Process. The upper loop adjusts the 
parameters of the regulator, in such a way that the error 
between the model output and the process output tends to zero. 
In this type of control system, the main problem is to deter- 
mine the adjustment mechanism so that a stable system results, 
which brings the tracking error to zero. 

Model reference adaptive systems can be further subdivided 
into two categories: direct and indirect. In indirect con- 
trol the plant parameters are estimated and the control param- 
eters are adjusted based on these estimates so that the overall 
plant transfer function matches that of the reference model. 
Mimaiirect Control no effort is made to identify the plant 
parameters but the control parameters are directly adjusted to 
Minimize the error between plant and model outputs. 

It turns out that the model reference approach is applica- 
ble only to plants with stable zeroes. In fact the only way 
the closed loop transfer function (Plant & Regulator) can 
match the one of the model in its poles and zeroes, is by re- 
moving the plant zeroes by cancellation with closed loop poles. 
This operation leads to the presence of uncontrollable or 
unobservable modes in the closed loop systems, which can be 
accepted only if they are stable (i.e., their effect decays 
fo zero with time), and this constitutes the major limitation 


for model reference adaptive systems. 


C. SELF-TUNING REGULATORS AND POLE PLACEMENT 
A third approach to the adaptive problem is the self- 
tuning regulator. A block diagram of this type control system 


is GiVEn In Figure 1.3. 
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Figure 1.3. Block Diagram of a Self-Tuning Regqulazer 


The reference model of the previous approach is replaced 
by some more general desired performances; such as error 
minimization or desired closed loop poles. 

There are two loops in the system configuration. The 


lower loop contains the plant and linear feedback regulator. 
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The upper loop consists of a recursive parameter estimator 
and a design calculation adjusts the parameters of the 
regulator. 

Also it is possible to classify this adaptive control 
scheme as direct and indirect. This classification depends 
on the complexity of the design calculation block that is 


seen in Figure 1.4. 
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Figure 1.4. General Block Diagram of an Adaptive 
Control System 


iol 


The first method is to parameterize the system directly 
in terms of compensator parameters. 

This research project will concentrate on the last method, 
that is, pole placement and indirect adaptive control. The 
evaluation of the control law is indirectly determined on 
the basis of parameter estimates. This method is also called 
explicit, since the design is based on an explicit estimation 
of the process model.. The parameters of the controller are 
updated indirectly via estimation of the process parameters. 
Several algorithms to estimate the parameters of a linear 
model are available in the literature. In this thesis two 
well-known estimation algorithms will be investigated: Recur- 
sive Least Squares and Projection. They represent a tradeoff 
between complexity of computation and performances, in the 
sense that best performances (with relatively high complexity 
are obtained by using Recursive Least Squares. Also, Block- 
processing is investigated to determine the period of the 
adaptation of the control parameters. Finite time persistency 
of excitation is also studied and simulated for indirect 
adaptive control. 

The main reason of choice of indirect adaptive control 
for this research comes from the fact that it is applicable 
to nonminimum-phase systems, and there are no limitations on 
zeroes of the plant. Therefore, it is more general than model 
reference adaptive control. 

Because of the effects of the zeroes of sampled data 


systems on adaptive control theory, we start investigating 


i 


the behavior of zeroes of sampled data systems. This thesis 
is organized as follows: In Chapter II zeroes of sampled 
data systems are analyzed by stating results available in 
the literature, and simulation studies. The indirect adaptive 
control problem is presented in Chapter III and computer 
simulation studies are given in Chapter IV. 

The computer programs are presented in the Appendices B, 
C and D. Also the description of the programs is given in 


Appendix A. 


if 


It. JZEROES OF SAMPEED DATA SYSTEMS 


Many control strategies are based on the assumption that 
the zeroes of the plant are on a stable region of the com- 
plex plane so that they can be cancelled by precompensation 
or closed loop poles. One example is the model reference 
adaptive control mentioned in the Introduction. However it 
turns out that in sampled data systems with Zero Order Hold 
(the most popular mean of Digital to Analog conversion) some 
zeroes of the resulting discrete time plant are often in the 
unstable region. In this chapter we analyze the position of 
the zeroes of sampled data systems, in relation with the con- 
tinuous time plant dynamics and sampling frequency. The 
structure, the number of zeroes outside the unit disc and the 
low frequency characteristic of the pulse transfer function | 
are examined from the transfer function for an important set 
of continuous time processes. 

Finally, it is investigated that zeroes of the sampled 
data systems are sensitive to high frequency poles present 
in continuous time transfer function. 

A. TIME DOMAIN ANALYSIS OF THE ZEROES OF SAMPLED DATA 

TRANSFER FUNCTIONS 

Poles and zeroes are important parameters of linear time- 
invariant systems. The zeroes describe the way the internal 


variables are coupled to the inputs and the outputs. As 


14 


described in [Ref. 6], the unstable zeroes limit the perfor- 
mance of control systems, since many design techniques are 
based on the cancellation of the process zeroes. This can be 
done provided the system's Zeroes are stable. 

The transformation of poles in the continuous time domain 


to the corresponding discrete time is given as: 
Doe eT (2) 


where T is the sampling period. This transformation maps 

the left-half part of the s-plane onto the unit disc, so that 
stability is preserved. But there is no simple transformation 
for. zeroes from continuous to discrete time domain. The type 
of hold circuit affects the position of the zeroes. Most 
digital control systems use a zero-order hold, and for this 
type of hold circuit, the effects are considered. 

In the following discussion, the main results are limit 
theorems, which give the zero locations for small and large 
sampling periods. 

If the continuous time transfer function G(S) is rational 
feetollows from Equation (2.2) that H(z) is also a rational 


monet ion 


H(z) = (1-2-7) [( (eit) /(z-eP47) Res, 
ae 


,G(s)/s (222) 


where: 
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Reo G(s)/s(s-P,) 
S-eu 

and P. are the poles of G(s)/s. 

The £uncrion H(z) nas geneereaiae, n-l zeroes. For particu- 
lar values of the sampling period some zeroes may, however, 
go to infinity, or they may be cancelled by poles, i.e., 
hidden modes. Hidden modes should not be considered as zeroes 
of H(z). In fact, there are in general no simple closed form 
expressions for the zeroes of H. The limiting cases for small 
Or large sampling periods can be characterized. That is ex- 
plained in the following theorem. The major steps in the 
proof are given by [Reto wie 

Theorem 1: 


Let G(s) be a rational function 


(s=z,) (S=z,) «1. (sez) 
Se eran its i an (a 
Py Do) --- Pn 
and H(z) the corresponding pulse transfer function. Assume 
that m <n. Then as the sampling period T + 0, m Zeroes tee 
H(z) go to l as exp(z.T) and the remaining nom] 7Zere= aes 


H(z) go to the zeroes of Bo-m'2) where BL, (2) is the polynomial 


defined as 
n + bon +... + b (2 


and 
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k 
se ee esi eee) ee les pn (2.5) 
C=1 


The following results can be observed from the previous 
theorem: 


1. The limiting zeroes of a pulse transfer function depend 
critically on the pole excess of the corresponding 
continuous time system. 


2. A continuous time system with a pole excess larger 
than two will always give a pulse transfer function 
with zeroes outside the unit disc provided that the 
sampling period is sufficiently short. This may 
happen for quite reasonable sampling periods, and 
sainpled data systems with unstable inverses are thus 
quite common. 


An example is given to illustrate the above discussion. 
Example 2.1: 


Consider a system with transfer function 


I. 


G(s) = —— 
fee) 


The corresponding pulse transfer function can be computed as 


b 2° + bi,z + b 


Rail 2 3 
a 8 ae 
(Z-e=) 
where: 
b, = 1- (+t ania e 
ay = aaa at p2jee 5 fey Sata oe 


1 


be = (ee Mie ice oe 


H(z) has a zero outside the unit disc if 0) <= 92) eee 

The result of Theorem 1 can be understood by studying the 
behavior of the continuous time transfer function G(s) =s ™ 
and its pulse transfer function. The reason why Theorem 1 
holds is because for a sampling period T small every system 


n-m 


behaves like G(s) = l/s with the number of poles n and 


number of zeroes m. 


The pulse transfer function corresponding to G(s) = a 
is given by 
7 BL {2) 
Bizgy -= aaa (236) 
Ts 2 ss) 


where B 2) is given in Equation (2.4). 
The polynomials Ba are found for some values of n using 


the program given in Appendix A. 


B, (2) — 
Bo (z) = z+1 
2 
B, (2) = 2 +4241 
ye = 2° ieee 
Sea 3° 4267° 266e> ean 


ILS 


The polynomials Ba HaveuazecOcsseutside Of on the unit 


wuaerce ror i 7 2 Themiunstable Zeroes are given in Table 2.1. 


TABLE Ziel 


UNSTABLE ZEROES OF Bf) 


n Unstable Zeroes of BZ) 
2 ~1 

8 =o s2 

4 ie He) ey 

B eee ee) 


Therefore, it can be noticed that there are continuous 
time systems with stable zeroes such that the corresponding 
pulse transfer function has unstable zeroes. 

Deets DOssitble to give a complete characterization of the 
zeroes of the pulse transfer function for small sampling 
periods. A similar result for large sampling periods is 
given by Theorem 2. 

Theorem 2: 

ie (s/s Oeudestricaly proper rational transfer function 
with G(0) #4 0 and ReP; < OS Then all zeroes of the pulse 
transfer function go to zero as the sampling period T goes 


emir inity. 


Ne) 


A lower limit on T for stable zeroes was obtained in 


[Ref. 14], here stated only for simplespolece 


Al 


Te = Qn[2A(n+1) ] (2575) 
where: 
6 = - max ReP. > O (23am) 
= aL 
and 
De S= max |A,/G(0) | (235g) 
7 
bf G(0) = O it follows from Equation (2502 
i P.?T 
= -1 -l e (2S 0e) 
H(z) = G(0)z + (1-2 ) es AS ERE 


The corresponding pulse transfer function has one zero at 
z= 1. The behavior of the rest of the zeroes may be more 
complex as is shown by Example 2.2. 

Example 2.2: 


Let continuous time transfer function be 


Ss 
Peete ee) 


Then the corresponding pulse transfer function is 
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He = e T(z-1){(e 1+ sinT -cos T)z aur ti ebwsin T+ cosT)]} 


Qz(z-e -") (27 Doe cosT eno, 
The zeroes are a = 1 and 
. = sree ee T+ cos T) et 
oe, -T 


e +sin T- cos T 


For control system design in discrete time domain it is 
important to know whether the zeroes of pulse transfer function 
are inside the unit disc or not. When all zeroes are inside 
the unit disc the sampled system has a stable inverse and all 
its zeroes can be cancelled. It is therefore of interest to 
find sufficient conditions which guarantee that all zeroes of 
a sampled transfer function are inside the unit disc. The 
criteria are given in Theorem 3 by [Ref. 6]. 

Baeerem 3: 


iiec(sj. lS a Strictly proper, rational transfer function 


with 
(1) ReP . — 0 
(rt) G(s) = 0 
Pees = arg Gliiw) < 0, for 0 < w < @ 


Then all the zeroes of the corresponding pulse transfer 
function H(z) are stable, so that criteria for no unstable 
momees dre given both in terms of G(s) and in terms of condi- 


tions on the Nyquist curve G(iw). 


Ele 


B. FREQUENCY DOMAIN ANALYSIS OF THE ZEROES OF SAMPLED 
DATA Sy oleh 


In this section, the number of zeroes outside the unit 
disc and the low frequency characteristic of the pulse trans- 
fer function are analyzed from the transfer function Of a 
set of continuous time processes. 

The discrete frequency response of a continuous process 
G*¥(w) can be derived from Equation (2.11) by substitute, 


. 


Z = exp(jwT) : 


(Z2-ceue > <2 — Cee) 
Gn =| es (21m 
(Z-p,) -++ (Z-p)) 


Or it can be expressed by the holding device transfer function 


H(w) and process transfer function G(w). Then it becomes: 
Gt(w) = = J H(wtkQ) G(wtkQ) (2 ey 
k=-00 
where k is an integer and Q = 2n/T. G(w) is the frequency 


response of the continuous time system G(s). 
For zero-order hold, the holding device frequency 
response is 
A yeres 


H(wt+tkQ) = 3 (wm +kQ) (2.230 


Substituting Equation (2.13) into Equation (2.12) and calcu 


into account that G(-w) = G(w), where G is the conjugated 
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value of G, we obtain 











de Bee 
G* (i) = 7 Sm Choy T= f(a) +e (a) } (Cee le) 
where: 
> Ww G(Q-w) 
1) = G=n (etn) C2 5)) 
_ ee Cire a) Celery) Gu). & 
ey acta) a on” (eee eee 


and r is an integer. 
The number of the zeroes outside the unit disc are given 
by Theorem 4 [Ref. 10]. 
Theorem 4: 
ifm cCONLinuous pr@eess satisfies Equation (2.16) and the 


@ondition 


G(2-w) 
the number of poles of G(z) outside the unit disc is zero and 
the phase angle $(Q2/2) = arc G(Q2/2) is between 0 and -n then 
its pulse transfer function possesses no zeroes outside the 

unit disc. If the value of phase angle is between -7 and -27 


one of the zeroes lies outside the unit disc, etc. 


ae 


(z-9)) cae (z-9,) 


Saye. T/T 
(Z-e ) ... (z-e 


eg) 
5 


) 

Thus, inverse stable continuous processes frequently 
possess inverse unstable pulse transfer functions. For 
instance if in a transfer function satisfying the condi tiem. 
given by Equations (2°16) “and (227 ae and T; > T are 
real and positive for.every i and the orders of its denominator 
and numerator differ by more than two, then usually thegpae 
transfer function is inversely unstable. 

Also, from Theorem 4, the phase of the continuous frequency 
response at the Nyquist frequency 2/2 is directly related to 
the number of unstable discrete time zeroes. High frequency 
dynamics (where high frequency is intended as above the Nyquist 
frequency) influence the phase at the Nyquist frequency, so 
that it can cause the number of unstable discrete time zeroes 
to increase. This high frequency pole may come from the struc=- 
ture of the measuring devices, actuators or other hardware parts 
of the physical realization. Usually, we neglect these terms 
in the transfer function. The following example can give an 
idea about this problem. 

Example 2.3: 


Let the transfer function of the plant be 





which has a pulse transfer function given by 
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Hez) = 


with one zero 


By adding 


slic 4 = 


D(s) = 


the zeroes of 


tem for tf = 


2 


woe U2 + 
2 


0.9 and some We 


amc plotted in Figure 2.1. 
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the extra dynamics D(s) as 
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values are given in Table 2.2 
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We notice that as a decreases, the unstable zero gets 


more and more unstable. 


ey 


III. INDIRECT ADAPTIVE CONTROL FOR DISCRETE Wi eoveoe = 


As mentioned in the Introduction, the main advantage of 
indirect adaptive control systems comes from the fact that it 
is applicable to discrete time systems with unstable zeroes. 


It is assumed that the system is as given in Figure 3.1, 





Figure 3.1. Single Input--Single Output Plant Sysrem 


The pulse transfer function is given by 


H(z) = (3255) 





r(z) and p(z) being polynomials given by: 


B4z) =) 255 Py2 2 ees (32928 


re) a nae (3539 


Ul 
a] 
ke 
N 
aa 
a) 
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The degree of the Jememina toe polynomial p(z) determines the 
system order n. Also, if we assume the plant to be causal, 
the polynomial r(z) has degree at most n-l. The sequences 
u(t) and y(t) denote the system input and output, respectively. 
An alternative way of describing the plant is by its differ- 


ence equation written in operator form, as 


jo OD AlG en Mince an 0) bale) (3.4) 

where: 
GN) hae 18 8) = a pd” (3.5) 
(Me, Fr (3.6) 


and D=q is the backward-shift operator, defined as 
ee) = y(t-1). 

In the adaptive control problem the parameters in p(D) and 
r(D) are assumed to be unknown. Only the system order n is 
assumed to be known to the designer. Hence, two problems 
appear at this point: 

1. Parameter estimation; and 
2. Compensator structure. 

The general block diagram of the indirect adaptive control 
problem is given in Figure 3.2. Estimation of the parameters 
of the plant is mentioned in Chapter III.C. The compensator 


structure is derived from the pole placement problem for 
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ESTIMATOR 





¥ (CE) 





CONTROLLER 


Figure 3.2. General Block Diagram of IndirecteAdapea 
Control 


linear systems; and it is computed on the basis of the esti- 
mated plant parameters. 

In the next section, we discuss the pole-placement 
problem which will be the basis for the determination of a 


Suitable controller structure. 


A. POLE PLACEMENT DESIGN BASED ON INPUT-OUTPUT MODELS 
The purpose of pole-placement by state-feedback is to 

determine a feedback controller so that all poles of the 

closed-loop system assume prescribed values. This can be 


easily achieved if all state variables of the plant are 
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measured. In general, the states may not be directly avail- 
able. However, under certain observability conditions the 
state can be observed on the basis of input-output measure- 
ments. For this reason observers are used to estimate the 
state of any observable realization. Along these lines, the 
controller can be considered as composed of an observer and 
eegain matrix. 
Consider a state-space representation of the plant assumed 


to be controllable and observable: 


x(t+l) = x(t) + Tu(t) aD 


y(t) = Cx(t) (3.8) 


with o. [; C matrices of dimensions n xn, n xl and 1 xn, 
respectively. 

The block diagram of the controlled system combined with 
meme OoSerVver—-Controller is Given in Figure 3.3. V{t) is an 
external input which will be discussed in Chapter IV. If 
we restrict ourselves to the class of linear control inputs 


we can write: 


u(t) = - Lx(t) + v(t) (229) 


as discussed in [Ref. 3], where L is a constant matrix as 


oul 





<j ESTIMATOR 






Figure 3.3. Block Diagram of the Controlled ss -cen 


and x(t) indicates the estimated values of the actual states 
of the plant. It is a well-known result in systems theory 


[Ref. 3] that we may define the observer dynamics as 
x(Gtl) -=> x(t) + Tule) + Kiy (2) exe (32 oF 


where K is a matrix such that 


and $-KC is a matrix with eigenvalues inside the unit circle. 
Rearranging Equation (3.10) we obtain the state-space 
dynamical equations of the controller as a two input-~-one 


output linear system. 


x(t+l) = [O-KC]x(t) + Tu(t) + Ky(t) (7p 


a2 


UD COM (3°12) 


To convert the state-space representation of the pole- 
placement problem to a transfer function (using polynomials) 
form, taking the z transform of Equations (3.11) and (3.12) 


we obtain 








X(z) = (21-0) tru(z) + (21-0) ~tKy(z) (mB) 
-1 =I 
Diz) = =L(ziI-O) TU(z) =L(zI-O) ~“Ky(z) +V(z) Ge 14) 
where we define the matrix Q = $-KC. Also we can define the 
rational functions 
7 eel » -Goaaq t2i-0)7 
Eel) NS det (z1I-Q) 
eee Cz } 
> ata Gee's) 
and 
= ae _ o-L adj (zI-Q)K 
ee ae det (z1-@) 
Es) 
aq(z) eis) 
where q(z), the observer polynomial, is an arbitrary stable 
polynomial. Arbitrariness of q(z) is guaranteed by the 


assumption of the plant being observable, while h(z) and 


S23 


k(z), the controller polynomials, have to be computed, on the 
basis of the plant dynamics and q(z)- 
Thus, the structure of the control input u can be described 


ass 


kz) 
oi) 


hz) 


a2) q(zZ) 


U(z) + 








Y (2) ee ee) (3 ae 
Or it may be expressed in a difference operator form as 
qiD)u(t) = Kip) ule e) ser ee) (32 18ee 


To make the notation more attractive, let us write 
Equation (3.18) as 


h(D) 


u(t) + aCe) 








yit) + ve) (3. 1o3 


The block diagram of the above controlled system is in Figure 
3.4. Hence, the closed-loop system defined from external 


input v(t) to the output y(t) has transfer fumeereg 


= GD) sD) 3 30 
eae p(D) Tq(b)-k(D) } =F (yA) V!*? (3 eee 
In the pole-placement problem, the goal is to determine 

the compensator parameters (k,h,q) so that the closed-loop 


poles are as we desire, 





y(t) = te v(t) (3m 
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Figure 3.4. Transfer Function Form of the Closed- 
Loop System 


p*(D) being an arbitrary stable polynomial in the difference 
Operator D. Roots of p*(D) are the new assigned poles of the 
closed-loop system. 

Bvecqudting  reuations (3.20) and (3.21), we obtain that 


k, h, q have to satisfy the polynomial equation 


Po igt(D) Seb) = ei(D)h(b) = q(D)p*(D) (3.22) 


which can be written as 


Rib J pip) aiD)e (D>) = FD) (oie) 


where: 


Et) eee — a CE) lo CD) = p* (D)] . 


a2 


Therefore, the problem of determining a suitable compen- 
sator for pole-placement is equivalent to solving the poly- : 
nomial equation (3.23) which is"known byeehe name OL 
Diophantine equation. The conditions by which this equation 
can be solved, together with the method of solution itself, 


are given in the following section. 


B. DIOPHANTINE EQUATION 
The general form of the Diophantine equation in the 


unknown polynomials K(z) and H(z) is given by: 


F(z) = biz pptz) = itz) ez (3.245 
with: 
m 
Liza Ko + kz to... t k 2 7 kn x 0 (3.225 
m 
hz) = h. + h,2 +... 4+ h2 (3.269 
ee n 
P(z) = Po t+ Py2 + «+ + DZ (32295) 
= n 
R(z) = ry + r42 +... t+ v2 (3220) 
and 
2 n+m 
Bez) = fy + £12 + £52 ob, ee erate iene (3.29) 
( 
where: 
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Operant a oer re and f. are constant, not necessarily 
all nonzero. 
In the general pole placement setting, the polynomials 
mez), F(z), F(z) are given, while h(z), k(z) are unknown. 
In the rest of this section, we determine the conditions 
by which the Diophantine equation can be solved and the method 
Gi solution. 


PmecubSstreuting Equations (3.25)—(3.29) into (3.24), we 


obtain: 
elec n 
i, +£,2 eatees tf ae = (Dy +p Zt... +p 2 ) (ky +k,z + 
m n 
+ k 2 ) + eam +r jz +... tr 2 yale + 
m 
h,z ae pees hZ ) (3230) 


and equating the coefficients of the same power of z yields 


the linear relation: 
a : l f., uae ie (3231) 
where the vector C is defined as: 


(Soe) 


ema the matrix San as: 


a7 


Fi + =:=poee 0 OU. 
P57 Pj eee O ee 
Pa re ee es 3c (333 
OG A'0 Po Py cee 
Oe ro rise. 





Equation (3.31) is a linear algebraic equation in the unknown 
vector C. The matrix San consists of m+l block rows; each 
block row has two rows and can be obtained by shifting its 
previous block row to the right by one column. It is a 
2(m+1) x (n+m+1l) matrix. 

In order for a solution to exist, the matrix oa has to 
be full column rank [Ref. 5]. This can be satisfied if 
2 (me) se nme or me oa. 

When =a is a square matrix (m = n-1), it is called the 
Sylvester matrix of P and R, which has nonzero determinant 
provided polynomials P and R are mutually coprime, i.e., they 
do not have common factors [Ref. 6]. 

We can summarize the steps to find h and k, mentioned 
above, as: 

1. Let n be the order of the system. 


2. Choose q(z) to be an arbitrary polynomial of dequeemm 
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cee corm the matrix Sin emo (Zz) and RUZ) < 


Zeeoce up tie polynomial F(Z) as: 
Ezy — oz) tp (2) = p*(z) ). 


PecOlVG £OL parameters of h(z) and k(z). 
Oe Set Up Ehe polynomials h(z) and k(z), or h(D) and 
k(D) using the relation D =z7l 
C. PARAMETER ESTIMATION 
As mentioned above, in indirect adaptive control, we 
have to estimate the parameters of the plant, where the 


parameters are the coefficients of the transfer function 


H(z) = (34) 





We can compute the controller parameters h(D) and k(D) 
from the- Diophantine equation on the basis of the estimated 
meant (say r,p) . 

A large number of different identification methods are 
available. In the literature, one broad distinction is 
Paeveen On-line methods and off-line methods. In the off- 
line case, it is presumed that all data are available prior 
to analysis. Consequently, the data may be treated as a 
Eemolece block Of infOrmation, with no strict time limit on 
mien process Of analysis. In contrast to the off-line case, 
the on-line case deals with sequential data, which requires 
the paraemter estimates to be recursively updated within the 


time limit imposed by the sampling period. 


As mentioned, on-line estimation schemes produce an 
updated parameter estimate within the time span between 
successive samples. 

Also, the on-line methods are the only alternative if 
the estimation is going to be used in an adaptive controller 
or if the process is time-varying. In this thesis, we con- 
Sider two classes of estimation algorithms: 

iy (projection; s 
2. recursive least-squares. 

Before proceeding, it can be said that input-output 

characteristics of a wide-class of linear and nonlinear 


deterministic dynamical systems can be described by a model 


expressed in the following form: 
ae 
y(t) =. o¢(b=1) oF (3.259) 
where y(t) denotes the system output at time t, @(t-1) 
denotes a vector given as: 
lk 
) 


d(t-1 = fy(t-l) y(t-2) ... y(t-n) u(t-1) ... u(t-n)] 


(3736) 


and @ denotes the parameters of the plant, described as: 


ae 


: ns {-Pyr Pore eer PyrYyr Tore - es ES (3238 


The following example illustrates the representation of 


the plant aS given ingequationgise 
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Example 3.1: 


Suppose the pulse transfer function is 


Bee) - z+l1 


Z +2zZ+4+3 


It is also expressed in the following form: 


- -2 
ne 1, e Z a2 = 
1+22 +32 
-1 -2 
1 tC ies: eee 
hig 7) = ~= Ses 


a ee sane 


Hence, the difference equation becomes 
ae ee Ee oy (Eo et) utt—2) 


which can be expressed as: 


pace — ly(t-1) y(t—-2) uf(t-1) u(t=2)] 





The rest of this section illustrates the estimation 


methods. The projection algorithm and the recursive least-squares 


algorithm are analyzed sequentially. 
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(1) Prejecticn 42lcenicam 
By the Projection Algorithm, the sequence of estimates 
8(t) is recursively computed as: 


n~n 


Ao eee ee ao (t-1) 


82 tt 2) ___[y (t) 0 (t-1) 78 (t-1)] (3.3188 
c +6(t-1)°6(t-1) 


with arbitrary initial estimate 8 (0) and c > 0; 0 < a =a 
The detailed derivation steps are given by [Ref. 2]. 

This algorithm is also known as the normalized least-mean- 
squares (NLMS) algorithm. The algorithm results from the 
following optimization problem: Given eet and y(t0, 


determine @(t) so that 


2 


i= File ct) -e(t-1) | | (3.39) 
is minimized subject to 
ape 
y(l)< (= oe aoc) (3.40) 


The main properties of the projection algorithm are the 
following: 
1. Estimation of 6 at time t, @(t) 1s alwaysieloescr ae. 


the actual value of 6 than the preceding estimated 
value-9o(f=1) . 


oe) oe a | | @ (t-1) =o < | | oC Scene (3.41) 
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Pee vwnen the time goes to infinity, the error between the 
actual output of the system and its predicted value 
will converge to zero. 


a ee sO (se) 


cae (he Gee eee pee 


Wwieinc: 
que 
e(t) = y(t) ~ ®(t-1)°6(t-1) 
Preeelaniiiie(t) —o(t—-K) || = © for any finite k (ds) 
{00 : 


Mie adaptation hate comverges to zero as £ > ©. 


ot) eee Least Squares Algorithm 

Becoceing Lo Gauss Che principle om which the least square 
estimates are based is that the unknown parameters of the 
Meeecess should be chosen in such a way that the sum of the 
squares of the differences between the actually observed and 
computed values multiplied by numbers that measure the degree 
Of precision is a minimum. 

The recursive least squares algorithm is given by the 


following equation in [Ref. 2]. 


‘i A‘ ee eet) 0 (eo 8 ead] 


cee = O(t-l) + a (3.44) 
ee Oe) Pe 2) O.CL— i ) 


mem tc > 1 and 
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BP (t=-2) 0 (t-1) 0 (eas een 


Gee (ee 7 
1 +O (t=) P(t 2 es 


(3.459) 


with 6(0) given and P(-l) is any positive definite matrix 
coe 
O 

The algorithm results from the minimizaticnee. men] 


following quadratic cost function: 


ut “ uw = - 
ar x (6-6 (0) Po (6= 9. Coe 


(ej ote 
1 


(Sa 
rap) 
lt 
No} 
teayg 


S 
(3.46) 


We can observe that the cost function represents the sum of 
Squares of the output prediction error e(t). 
Pa aus) ses ise) ys (3.47) 

The second term of the right-hand side of the cost function 
accounts for the initial parameter estimates weighted by the 
matrix Bare 2 this way, we can consider Ee (the initial ome 
dition of P(t) in Equation (3.45)) as the “conttdene aan 
the initial conditions 6(0). 

The least squares algorithm, as can be seen from the 


Simulation results, has much faster convergence than the 


ProOyection algorichin. 


D> PERSISTENCY OF EBXGCETALTICH 
In order to guarantee global stability of indirect adap- 


tive control algorithms, the input to the plan= hace ome. 
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persistently exciting, which in turn implies global conver- 
gence of the plant parameters to the corresponding actual 
values [Ref. oe To obtain consistent estimates of the plant 
parameters, it is necessary that the input signal to the plant 
be sufficiently rich in frequencies, and excites all modes 
@eethe plant. 

In a closed loop set-up, the control input is the sum of 
an external input signal v(t) and a feedback signal from 
the adaptive controller. 

The feedback signal may in principle cancel any excitation 
contained in the external input signal. This problem and 
the potential for unbounded growth of the control and output 
Signals have made the guarantee of persistent excitation a 
G@ifficult problem. Recently, Elliot [Ref. 8] gave sufficient 
conditions which guarantee persistency of excitation. The 
following theorems summarize these conditions. 

Theorem 4.1: 

Let w be the number of parameters estimated. In order to 
guarantee global convergence of the plant parameters to their 
true values, the following conditions have to be satisfied: 


1. The external input v(t) should consist of a sum of 2w 
sinusoids. 


2. The compensator parameters should be updated each N 
samples; with N > lOn- 


Therefore, in this thesis report, block processing is 
used in the sense that N data samples are taken, and N iter- 
ations of the recursive least-squares algorithm are performed 


nN NAN 


between control parameters (h,k) updates. To avoid time 
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variation of the plant dynamics during the spanning process, 
the control parameters are held constant during spanning 
blocks (of length N) and are changed only between them. 

As we can see from Theorem 4.1, these two conditions 
are sufficient to guarantee parameter convergence for any 


ine tal cCondie toms ase 


ati | — eG 
kro k 


Oy being the estimate of 6@* at time ty: 
In the next section, we will examine finite time persis- 


tency of excitation. 


BE. FINITE TIME PERSISUENCY OF VEXeCr ration 

Tt is clear that the persistency of excitation condition 
on the external input v(t) is the main limitation about this 
adaptive control algorithm. Recently, in a report by Cram 
[Ref. 15], a possible solution to this problem has been given, 
by stopping parameter adaptation when the performances are 
close to the desired ones. More precisely, the model output 
error between desired output of the closed-loop system and 
controlled plant's output is the measure of how far the system 
is from the desired performance (i.e., pole placement) and 
adaptation is stopped whenever the error falls below an arbi- 
trary, preassigned threshold. The configuration of the system 


is Given in Figure 3.5, 
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SWITCH DECISIO 






Figure 3.5. Indirect Adaptive Control Scheme with 
Added Finite Time P.E. 


The model output error can be defined as 


r, (D) 
Mees = yit) = sem (3243) 


The finite time persistency of excitation algorithm can 
be given as follows: 
me tc = ttl. 
mee compute y(t) from Equation (3.48). 


3. If {u(t)| > €, compute compensator parameters as in 
eiiapter il E. 


Bee be ju(e) | < € go to 1. 
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A major difficulty with this algorithm is to bevableuce 
guarantee global stability of the whole system. Therefore 
we have to prove that the signals in the loop are bounded, 
provided the external input v(t) is bounded. Global stability 
of the closed-loop system is proved below. 

When time goes to infinity and persistency of excitation 
is present, parameters of the estimates of the plant poly- 


nomials r, and P,_ converge to the actual values r and p. 


c 
Therefore, the difference between measured output and desired 
output tends to zero, and also H(t) tends to zero. Thus, by 


definition of limit, there exists an instant ty such titan 


PGs) || Be 
for all t > oe and for any given ¢€ > 0. Rearranging Equation 
(3.48) yields 
r(D) 
= + 3.49 
y (t) Say ee meee) ( ) 


which can be expressed on a block diagram form as in Figure 3.6. 





Figure 3.6. Block Diagram Representation of Equation 
(a42) 
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INE @eyo [6 cieeigepesinsiel ei ee input single output 
linear stable process with bounded input v(t) and bounded 
eaiaeuboance w(t). From the fact that stable linear systems 
with a bounded input produce a bounded output, the global 
stability results follow easily. So that with finite time 
persistency excitation, the closed-loop system is eventually 
equivalent to a system with poles as desired, bounded zeroes 
fete aA bounded output, disturbance u(t). 


We can generate the external input v(t) as 
vlc) = Vv, (t) + vate) (2250) 


with V(t) the external desired command and ed an added 
persistency excitation signal. When the adaptation is con- 
tinuing, the external input is composed of Ve Ce) and VEE), 
but after stopping the adaptation persistency of excitation 
is not needed, it will be identical to the longer, so that 
Vv, (t) can decay to zero. 

iieents way, the adaptive algorithm is activated only 
when the performance index u(t) is larger than a minimum 
threshold. In the next chapter, simulation studies will show 


their efficiency via some examples. 
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iV. SiMULATION VS LuUDmes 


This research report includes three computer programs. 
The program named CONDIS, given in Appendix B, is used to 
investigate the behavior of the zero in sampled data systems. 
This program computes the parameters of the discrete time 
transfer function from the parameters of its continuous time, 
and the sampling rate. This program has been used to inves- 
tigate the perturbations of the zeroes for different values 
of the sampling interval. 

The other two programs given in Appendix C and D simulate 
the indirect adaptive control using recursive least squares 
and projection algorithms. Entering order and numerator, 
denominator parameters of the aren observer polynomial and 
desired closed-loop characteristic polynomial, the program 
simulates the entire system in an interactive fashion. Also 
it is capable of graphical and tabulation results. 

The adaptive controller presented above has been simulated 
for several different plant dynamics. 

Example 4.1: 


Let the discrete time transfer function of the plant be: 


H(z) = 22 2 eee 
oe 27 +0.75 (2 =1. 572 = 075) 
The plant has one unstable zero, z = -2, and two poles 


P) = IE 3) eno fe Po = 0.3. 


210) 


Hiecmpelynomltalseq (zy and p*(z) are chosen to be stable 
Bolynomials with degree n. 


Pieper rienlanr,. let 


aa = 2° See t= Ore Bene tz =O. 7) (2 +0. 4) 


and 
pee — Le OSs = fz = 065) (2 ~- 0.6) 


bem Having stable roots, i1.e€., inside the unit disc in the 
Z-plane. Before going into simulation study, the problem is 
solved analytically by comparing their responses. By writing 


the plant dynamics in shift-operator form, 


regante 
-] a gq 


il =De - 30. een 


we can derive the difference equation of the given system as 


eee eae Oe oy (t-2) + W(t—1l) + 2u(t—2) 
y(t) = @(t-1)T9* 
where 


onl 


-y (t-2) 
@(t-1) = 
u(t-1) 
u(t-2) 
and 
22 
OS75 
co = 
: i 
2 


8* being the vector of the actual parameters of the plant. 


Generally, it can be partitioned as: 


"Oo 


p and r being vectors of parameters of the denominator sana 
numerator of the plant, respectively. 
From Equation (3.23), considering k{z) -and@int2) @a- 


controller polynomials, we can write 


In this example we have chosen: 


5(z) = 922°] oes 


2 


rz) = 2s 2 


GZ) a oe. 37 = O28 


Se2) eae = Vet OS 


By the constraints on the polynomials for solvability of 
the Diophantine equation, we choose k and h to be first order 


polynomials as: 


[CZ = kj2 + Ky 


Te) 


h,zZ + hy 


Substituting these polynomials into the Diophantine 


equation in matrix form yields: 


Ko t Ge << ets Bees | 
h, 2 ioe” °O O11 
ky ins =o | O72 

| | 
h, 0 ae Sn eG -0.9 i: 
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Since the plan does not have any pole-zero cancellation, 


the determinant of the matrix Sin is not equal to zero. Hence, 
it is solvable. Solving the above matrix equation, we compute 
the control parameters ki = -.899, Ky = 07 occ hy = —0 sca 


and hy = 0.195. The corresponding polynomials of the con- 


troller are: 


k (z) =. 0972, Ono 


hig) 0.39z + 0.195 


This corresponds to the optimal choice of compensator 
parameters for the desired pole placement. Therefore the 


compensator given by the difference equation 


co) = SO ea 


yields closed-loop transfer function 





ta) = 3g 
= (z)r (Zz) 
p(z) (q(2) -~haSe eie 
which becomes: 
H(z) = 24. 
Z =1.d2 40.3 
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The step response of this system is given in Figure 4.1. 
The next approach to the problem is simulation of the system 
using a computer program. It is observed that from the simu- 
lations, after the 10n iterations, the estimated parameter 


~ 


will be closer to 98* and 9 will tend to zero, where 


egies e= =| 


nm “nw 


This means that p and r are converging to the actual param- 
eters, and 7 and K converge to the actual values for the 
desired pole-placement. The desired output and controlled 
system's output are given in Figure 4.2. Also, parameter 
error and output prediction error are given in Figures 4.3 
and 4.4, respectively. After convergence of the prediction 
error to zero, the persistency of excitation added to the 
external command is turned off. In this example, we have 
chosen the step as external command. External input and plant 
input are in Figures 4.5 and 4.6, respectively. 

If we compare Figures 4.1 and 4.2, analytical and simu- 
lated system's results are almost the same. 

In the next example, it is considered that one of the 
plant parameter is changed after some period of time. It is 
investigated how the system behavior will change in this 
Sondition. 

Example 4.2: 
Let the plant be the same as the one in the previous 


example. After two blocks length, the zero of the plant is 
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perturbed from z = -2 to z = -3. To compute the parameters 


associated to the perturbed plant with pulse transfer function 


Z+ 3 
2 
Vi se Pee Ss (Chay 5 


write its difference equation as 
y(t) = 2y(t-1) - 0.75y(t-1) + u(t-1) + 3u(t-2) 


and immediately 


By using the same observer and desired polynomials g(z), 
p*(z) and following the same steps as in Example 4.1, the 
controller parameters become Ky = -0.894, Ko = -0.778, 


hh. = -0.305 and h. = 0.152. 
aL oO 


Then, the polynomials k(z) and h(z) are: 
K({z) = -=2.8942)- 107776 
h(z) = -0.305z + 0.152 


The closed loop transfer function becomes: 


4 


H(z) = a rr 
gay Dee IN ae Oa 


The step response of this transfer function is given in 
Bagure 4 Ps Also, the other graphical results are in Figures 
4.8-4.12. 

In the next example, we investigate how the disturbance 
at the output will affect the behavior of the controlled 
system. 

Example 4.3: 


Consider a plant with pulse transfer function as: 


H(z) = ———— 

Ca 27 tO. 75 
Pae@massSume an Output disturbance exists. Let the disturbance 
be sinusoidal with frequency 51/2 rad/sec. In this case, it 


is observed that the estimated parameters of the plant don't 
converge to the actual parameters. Corresponding plots are 
given in Figure 4.13 through 4.17. However, if the disturbance 
is small, the parameters converge close to the actual values 
and we can still obtain satisfactory performances. 
Comparison of Different Estimation Techniques: RLS and P.A. 
In the above examples we used a recursive least-squares 
algorithm for estimating the plant parameters. In this sec- 
tion, we compare the behavior of the indirect adaptive control, 
Maem the projection algorithm is used. 
The projection algorithm has the advantage of requiring 


less computations. Approximately, the number of computations 
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grows as ne in the recursive least-squares algorithm. But 
oaethne projection algorithm, it grows as n. As noted before, 


the projection algorithm can be described as 


a&(t-1) 


0 tern ly (t) -o(t-1) 76 (t-1)] 
Gap esi) a eae) 


e(t) = 6(t-1) + 
with c > 0 and 0 < a < 2. The value of the parameter a is 
crucial. The behavior os the algorithm greatly depends on 
mes value. 

In the next example, the projection algorithm is used 
for estimation procedure on the previous process. 

Example 4.4: 

Using the same pulse transfer function and controller 
polynomials as in Example 4.3 and choosing the constants 
a = 1 andc= 1 in the above equation, results can be observed 
from Figure 4.18 through 4.22. 

Actually, using this algorithm, several systems have 
been simulated. This result is the most reasonable one. 

Considering the above results, indirect adaptive control 
is a very effective control technique when the parameters of 
the plant are unknown and large and unpredictable variations 
occur in the plant dynamics. Also, it is applicable to 
nonminimum-phase systems. This can be considered as a big 


advantage of the method. 
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APPENDIX A 


DESCRUPULONSOr COMPUTER PROGRAMS 


Computer programs are written using WATFIV and FORTRAN 
programming languages. All of them are prepared as an 
interactive program. The following explanation is about how 
one can use these programs. 

The program named CONDIS given in Appendix B finds the 
pulse transfer function from the continuous time transfer 
function. The program asks for the sampling interval, orders 
and coefficients of the continuous time transfer function 
polynomials. The continuous time system is described in state 


variable form as 


and the corresponding discrete time system as 


kt) = ox (k) + Tu(k) 


This program gives the 6 and [ matrices corresponding to A 
ed B. 

There is no limitation about the order of the system in 
all programs. For higher order systems, the dimensions of 
the declared matrices should be increased. 

The program named RCIOP given in Appendix D simulates 


an indirect adaptive control system using the projection 
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algorithm. The program asks some questions to the user in 
the following order: 

l. The order of a plant (N); 

2. The constants a and c that are given in Equation (ym 


3. Coefficients of the numerator polynomial of the plant 
(in ascending order of 2); 


4. Coefficients of the denominator polynial of the 
plant (in ascending order of ez 


5. Coefficients of the controller polynomraleg (Zia 
D* (2s 


The external input is defined inside the program. It 
can be changed when it is necessary. 

The last program RCIOR given in Appendix C is the most 
important one. Least-squares algorithms are used for parameter 
estimation. This program asks the same questions about system, 
except question #2. 

Computer programs RCIOR and RCIOP give the graphical 
and tabulated results. Graphical outputs are obtained using 


DISSPLA. 
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APPENDIX B 


COMPUTER PROGRAM CONDIS 


WATFIV ONUK, XREF, EXT 
NEECER NM, oa HOES 
INTEGER N Piezo, Jy oni MZ NG 


See 
REAL Bh ae ID 2 ES On 10 
il 
7 
7, TRB 


8h 8 
ae OMAIN PULSE TRANSFER' 
,, EYNCTION WHEN GIVEN THE S-DOMAIN TRANSFER FUNC. 


PRINT, 'THE S-DOMAIN TRANSFER FUNCTION SHOULD BE 
ihe tHts FORM" 


H(S) <NUM(S) /DEN( S 1 : 
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16) At? TAR 


NUM 8 is T(MYSE ENS Lt sees TZ 1st 2) 
oe +D(N)S**N-1+. 21 24+D(2)S+D( 1 
PRINT, 'ENTER THE NUMERATOR POLYNOMIAL DEGREE' 
BRING, T 


BEAD ENTER THE DENOMINATOR POLYNOMIAL DEGREE' 


M2=Mi+1 
DO 2 11=1, 
PRINT at ta yee” 
READ, TT( 115 
CONTINUE 
N3=N-1 
M3=M1+2 


TR (Mi Lt NS) THEN 
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12 


22 
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bop 
PUI 


DO 82 I11=M3,N 
Trl) =o a0 
CONTINUE 
ND IF 
DO 12 11=1,N 
PRINT ae ka 11, j=) 
READ 


coneene ee 
He 1i=1,N3 


=) 


/N 
56 Yee: T1=1,N3 
IF (J.EQ.13) THEN 
Aor: 
AA( 11, J)=0.0 
END IF 
CONTINUE 


CONTINUE 
CONTINUE 
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CALL FAMILY 
CALL RESULT 


tae sna) Wil ia) tWid2, 10)) 


TH4= “thi /i6 
PRINT, T 
IF (TRI. LE.TC) THEN 
TR2=TR1/2 


R1/9 
TR1O=TR1/10 
TRII=IRI/11 


CALL SCALAR 
CALL CALCUL 
CALL SCALAR 
CALL CALCUL 
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CALL CALCUL 
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CALE Gancu (N, NUN, 1, kh Dilwaay aay) 
CALL SCALAR( NN, 1, kK, TRo, Disco) 
TR4=( TR**4 
CALL CAnCU 
CALL SCALAR 
TR5=( TR**S 
GALE CAncU 


K, D2, AR D1) 
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Kp 
TR 
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IASB 310 
CALL CALCUE(N,M,N,1,K Ro EE eee 
CALL CALCUL se M,N, 1K, kK, GG7ee ya 


Rint fete 
PRIN 
PRINT, ' NUMERATOR COEFF.S OF THE PULSE TRANSFER 


, EYNCTION 
PRINT, 

PRINT, 'NUMCOEF' ,CO(N) 
BO 14700190 

EXECUTE AR’ 


J,PP,FI,RR) 
LP) ID,Ss} 

S RR) 

_J,EE BeoGn 
1, GG) @@ern 


L=} ,N3 

PRINT, 'NUMCOEF' ,CO(L) 
CONTINUE 

PRONG o 

PRINT 


PRINT, 'DENOMINATOR COEFF.S OF THE PULSE TRANSFER 
FUNCTION 


PRINT,' 
DO 151 L=1,N 
PRINT, DENUMCO' , ALPHA( L) 


STOP 
REMOTE BLOCK AR 
GAnn CALCUL(N, N,N,1,K,K, PP, FI,RR) 
DO 2480 1—1) 

DO 149 5s i, a 

mee efaR EN 
ARE ARIPEEL ab) 
CONTINUE, 


CONT INUE 
Be BLOCK 


EN 

***THIS SUBROUTINE READS THE MATRICES FROM THE DATA 
SUBROUTINE ENTER( J, D,/ EAE Ve) 

NTEGER J EB’, EF 


TE 
READ ClaoD 


’ é ‘ 
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pO READS ®S0/ Ic E, F), F=1,D) 
Eerie Sa 
ONT INU 


oul 
ele) 


END 
e ***THIS S SUBROUTINE PRINTS THE ARRAYS AS AN MATRIX — 

SUBROUTINE oCQPE(H, ©, 12,8, Gil, ea) 

INTEGER 

REAL bia ON iC, 0) 


Dd 132 $1.0 =K1(P,S) 
CONTINUE 

CONTINUE 

RETURN 


END 
SUBROUTINE RESULT(H, O,P, 5,7) 
INTEGER H, 0, P, S$ 

REAL T ( H, O 


PRINT 70,(T cen) 

poan( 6!" CB 74k 10( 8x, F12. 6)) 
CONTINU 

RETURN 


END 
C ***THIS SUBROUTINE CALCULATES THE MULTIPLICATION OF 
TWO MATRICES 
SUBROUTINE “cALCUL( u V,Y,Z,X,ZX,2ZT,W,Q) 
INTEGER U,V 
REAL 2z(U, Vv), wwe *oh Q(Y,V) 
DO 93 ou 
DO 792 RE1.,V WV 9. , 
“DO. giz Y 
T(2) x ee Zor Z, X)+W(Z,ZX)*O( ZX,X) 


bf 
WwW 
Hh) 


O~] 
OO 


el cont 

a2 CONTINUE 

93 CONTINUE 
RETURN 


END 

C ***THIS SUBROUTINE, WRITES THE MATRIX NAME AND ROW, 
COLUMN NUMBER 

SUBROUTINE BAMIT. ¥( ROW, COLUMN, NAME ) 

INTEGER ROW, COLU 

CHARACTER* 2 NAHE, 

PRINT 94 MATRIX 
94 ee 


ion Fa, A6,T19,1X,T20,A2) 
PRINT 95, ‘ROW NGuBES OW 


89 





95 
96 


QQ 


EORMAT(T1, No ee 10x ,T12,Al11,123,8X,T31, 12) 
PRINT ' ¢OROMN NUMBER. U 

FORMAT( TA, P21 Ti, ee Me SX ,13 1, 12) 
RETURN 


END 
***kTHIS SUBROUTINE CALCULATES THE SUMMATION OF THE 
absolute values OF THE SAME SOLAN EEEMENiS. *** 


SUBROUTINE ALI 1 12,H1,H2,G1, 
INTEGER 11,12 2 
REAL W2(1, isu, a 


42 H1=1,11 
W2¢ 1, HD) =W2( 1, H2 )+ABS(G1(H1,H2) ) 
CONTINUE 
CONT INUE 
RETURN 


END 
Cc eeu eR UTINE MULTIPLIES THE MATRIX BY SCALAR 


44 
43 


SUBROUTINE SCALAR( KI, Koy Glico Pics G4) 
INTEGER K1,K2 
REAL G3( Kl, joe, K2),P1 


bo. 44 Gil 
G4(G11, G2) =63(c11, G2)*P1 
CONTINUE 
CONTINUE 
RETURN 


END 
Ce ** THs UE tr CALCULATES THE SUM OF THE TWO 


TR 
SUBROUTINE put ks k K4,15,K5,X1,X2,X3) 
INTEGER K3 
REAL K1(K3 Ra atk, K4) ,X3(K3,K4) 


DO 46°15 
X3( 15, Ke )aKICIS, K5 )+X2(15,K5) 


46 CONTINUE 
45 CONTINUE 
RETURN 
END 
SENTRY 


90 


+ OF 


* 


QAQAQAAAAAANAANQ 


Ty 


8 
20 


APPENDIX C 


COMPUTER PROGRAM RCIOR 


KRHEKKKKKKKEKKKKRKAKRKRKKRKRKKKKKKKRKRKRKRKRKRKKKKKKKKEK 


THESIS PROGRAM 

INDIRECT ADAPTIVE ov 
PARAMETER ESTIMATI 

USING REC. LEAST SQUARE ALGORITHM 


+t + tt 


KEKEKKKKKKEKKEKKKKKRKEKEKKKKKKKKKKKKKKKKKKKKKKKE 


RRKEKKKKKKEKEKKKKRKRKEKRKKKKKKRKKKKKRKKKKRKKKKKRKRKRKEEE 


DEFINITION OF VARIABLES: es 


KREKKKKEKKRKAKKKKKKHKKRKRKRKRKRKKKRKKRRKKKKKRKKRKKRKRKKREK 


INTEGER N1,N,L,K,L1,1,M,02,30,No,N2, 41) seo 
INTEGER N4/N5_N6. COLUM, Js, m5/J57,J38, J3, L5,M6,L7 
INTEGER Bat 2 bEcs, ques ‘M7,L8.M2 Bey Pet M8 


vats 
reat valet rah aa ta iti 
REAL P4 


4 Ht ee gon 50} 


1903 #36163) 


Seccatae 


CONTINUE 


M= 
KRERKEEEKRKKKKKKKKEKEKRRKKKKRKKRRRKKEKERKERAKKKKRKKESE 


6: 
C * THIS PART OF THE PROGRAM ASKS THE sis late 
C * ORDER , PLANT AND MODEL NUMERATOR, DENOMI. * 


ell 


C * POLYNOMIAL COEFFICIENTS. ~ 
Cc KRKEKKKEKE 


RRREKKEKREEKKEKKERER 


7p 


2002 


2003 


2004 


Z001 


2005 


2006 


2007 


2009 
2008 


2 


243 


124 


“FORMA 


+ 


RREKKKKKKEKKKKEKKKEKEKE 


WRITE( 5,701) 
({1"q ENTER THE ORDER OF THE SYSTEM") 


M7=10*N*2 
L8=10*N2 
CALL FRICMS( 'CLRSCRN ') 


WRITE( 6,200 
FORMAT ( NTER PLANT NUMERATOR POLYNOMIAL 
COEFFICIENTS") 


WRITE( 6,290 
FORMA N ASCENDING ORDER OF 2Z') 
DO 200 I=1,N 


J7=1- 
WRITE( 6, 2004 


¢ d I 

FORMAT ( OEEE OF Z2<=( 2 .\)=—" ) 
READ( 5, * ) €5(T) 

CONTINU 


CALL FRTCMS('CLRSCRN ') 


WRITE( 6, 200 

FORMAT ( NTER_ PLANT DENOMINATOR POLYNOMIAL 

COEFFICIENTS" ) 

WRITE( 6, 2006 
N ASCENDING ORDER OF Z') 


WRITE(6, 200 
FORMAT( rol IGHEST COEFF. SHOULD BE 1. ') 
DO 200 I=1,N 


J8=I- 
WRITE 6 , 2909) J8 
FORMA @HEF. OF Z**(',12,')=') 
READ(5, *) C6(T) 

CONTINU 

Domi en= 


e260 
VP(L)=SIN LOB) SOME Ss UN ee 
+SIN( L*P1I/7)+SIN( L*PI/2)+SIN( L*PI/6 


CONTINU 
Bel 2431) 200 
VC(L)=100. 0 
CONTINU 
DOr i242 9r-1 
vin) —o0 0 
U( L)=VP(L) 
CONTINUE 
¥(NI =9. 0 
D TO 1=i 


a2 


CONT NUR Neer ee ae 
DO aos 1=1,N 


2010 
U 
2011 CONTINUE 
DO 


i CON 


eeniz 
TPC 1se a THEN 
P K)=30. 0 
BOC T,K)=0. 0 
nD {F ; 


CONTINUE 

CALL FAMILY(N2,N2, 2"), 

CALL RESULT (NZ, aN2 ‘ Pk 2. 

CALL COPYM( Ri ABs 

DO 123 I= ‘. Nee 

Ol( tM) —Onc 

oe CONTINUE, N20 le nn) 

s FERSERK im 


QQHH 


CAEL ERTCM NS ( 
702 EORMar( hy vane Q(Z) IN DESCENDING ORDER OF 2') 
703 ee ) HIGHEST COEFF. MUST BE 1."') 
704 Banoe (7?) te ee SHOULD BE STABLE POLYNOMIAL' ) 
DO 157 1-1, N 
WRITE( 6, 705 )N- ee: 
705 FORMA Ct COREE. OF Zee” i) 
READ(5, re) C3 
157 CONTINU 


CALE FRIGMS( "CERSCRN ") 
So 6 
FORMA 


706 fi ENTER PSTAR(Z) IN ASCENDING ORDER 
WRITE 

707 nies arcs | BIGHEST COEFF. MUST BE 1.') 

708 AMAT! EME ree IS THE DENOMINATOR 


28) 


709 


710 


7alal 


Neys) 


241 


ye) 
U107 


DO 


POLYNOMIAL' ) 
WRITE( 5,709 


Cn log MODEL ' ) 
WRITE(S , 7 
@ a -2:PSTAR(Z) SHOULD BE STABLE 


soyncatt 


al 

BITE B' 4 1 

FORM E(6, t COEFF. OF 2*4(77 12, )=') 
MBAS 8. pe) C40 1) 


114 J=N1 
CALL CNLOUE N2,M, N2 ike hoe Po eT 
GALL CALCUL( N MoN2, 1,L,k, V2, 25) iy? 
CAEPMCALC aM! Men2, 1,b,K, Vinee ul) 
sedate -V4(M,M§ 
Site SCALAR(N2 ,1 it 
A Sane te 
CALE CALCUL MM N te zs ne Kea , FT, £3) 
TSM=1. /( 1+P (ff 
CALL aA CUL io 5, ET, B2 
CALL CALCUL( N2, Nt Re 
N2, a! , Ly P7 FDSe 2 
Nz i, fsM bs) 
CALL SUB N2,N ae sai b2' Ps, 


K, TSN, V1,V5) 


2a) 


i 


,P 
P6 
st 67 





IF(J.EO.10) THEN 
CALL COPYM(N2,N2, Tyke, P2,AB) 


ean cebu Re st I ae P2,AB) 

CALL eG sees N I AB. NP2, AB) 
1e(J THEN 

cane COPYM( N2, id ae a2. AB) 

Cann Hae N i “ite NP2, AB) 

CALL Coby R2° oe iy Ee NP2, AB) 


So ahee NZ) Lp L PZ eo) 


94 


148 


241 


212 


74,9) 0Ual 


4001 


7007 


6503 


ELSE IF EQ, 7 
CALL CO vin ve ,P2,AB) 


see LE.M7) GO To 9591 


N34 PEN 


U1(1,M)=0. 0 
ConTINUE 
CALL COPYM(N2,M,1,[,U},U2) 
CALL FAMILY(N2,N2, "62 
CALL RESULT( N2,N2/1,K,P2) 
DO 148 I= =1,N 


DO 271 I= 


=U1(1,M 
CONTINUE (1 ) 
DO 212 I=1,N 


B(T)2U1( L4, M) 

CONTINUE 

DO 7011 I=1,N 

O11) Le es 
Thue 


N 
1(N1)=1. 0 
2(N1)=0. 0 
> O3Cr tt) -C4(1) 
contaNun 
DO 7007s 1—i0t) 

03( 1) =03( M8) 
cont sNot 
DEG1=N+1 
DEG2=2*N 
DEG3=DEG2+ 
DO 6001 I= 


DEG2 
i 6020 0), OR. ((1=K)36n. DECI 


Ro¢ }, G0, 20, 6503 
R( CT) SR 12 Shey (ket) #2 (I-(K-1)) 


hs 


6002 CONTINUE 
6001 CONTINUE soot). 
Do $003. 121° NG 


BS R( K9) 
9003 CONT E 
N4=2*N+1 
M1=N+1 
DO 41 J4=1,M1 
DO 51 _ I=1,N4 
ea 4 
Aly J 


ray ee or oe (K.LT.M1)) THEN 


PUI 
HH 


ont 
CONTINUE 
DO 61 J4=1,N 
De Val eal 
ACI/N pete Sa)= =F 


IF 
me "GP, 0.) AND Sek 
CONTINUE 
CONTINUE 
DO 10 COLUM=1,N4 
SUM2=0. 
DO 70 J4=COLUM 
SUM2= =sun2+( A 34, COLUM) ) **2 


CONTINU 
SS ASICHA 
A= 


(K. LT.M1)) THEN 


ox! 
a 


=hB9( SicMa MA} 
LT.0.1) GO TO 6666 


ie TEWRITECE. 6 
C14 eae FORMAT( i *) s waTRIX A IS SINGULAR’ ) 
DO 80 nae = COLUM, N4 
El Ua ere) 

80 CONTINUE 

Siiigede) ot 

SUM2=0 

BO 9Osu4=coL, 

SUM2= SUMO *( AC J4, COLUM) -G( J4) )**2 


90 CONTINUE 
DO 100 J4=COLUM,N4 


96 


lo 
Or 
OO 


a0 


6666 


DO 110 K=COLUM 
IF LEME Bs(SuM2). LT. 0.00001) THEN 


Fie ; ote -G( J4 *( A( K, COLUM 


BP ae ts ve THEN 
EMP=1. +TEMP 


i 
H( J4, K)=TEMP 

CONTINUE 

CONTINUE 

Do 131, J4=COLUM NA 


OLUM, N 
TEMPSO” 
DO 140 L=COLUM 
TEMP= TEMPSH( 34, L)*A(L,K) 
CONTINUE 
c( 34 g4 ak) =TEMP 
CONT 


CONTINUE, 
DO_150 J4=COLUM,N4 
TEMP=0. 


DO 160 L=COLUM 
TEMP= TEMP SHC 4, L)*F(L) 
CONTINUE 
se J4) “TEMP 
CONTINUE 
DO 170 J4=COLUM, N4 
DO 180 K=COLUM, N4 
CicWe eel buliie 
CONT NUE 
Ri j4=00 74} 
CONTINUE 
XCN4)SEUN4) /A(N4 N4) 
I=N4-1 
SUM2=0. 
M2=I+1 


DO 30 J4=M2 
nous =SUti2+A( 1, J4)*X( J4) 


6 ie UryymEt 3! 66601) THEN 
END te 


=[= 
LET. lesa, i} GO TO Gia 
Ea) =e 


ay 


141 


142 


5001 


3003 


eo 


2508 


3004 


S005 


120 


Ue at Ul 
I] Ou Ob 
Z 
ram 


tiBC a)= er eSB CST see Ta at i sa(1)*y( J-1) 
Cees J tVeI)) 7 i=s201)) 


Do” 3903: ts 1,N 


cont tng 

DO 189 T=1,N1 
YM( 1)=0/0 

CONTINU 

YM( 35 )=0. 

DO’ 2508 I=1,N 
B1( 1) =c5( M8) 

conPINuE 


ac IMC IS) =YMC IS ) - -C4( 1) *YM( J5-1)+B1(1)*V( J5-1 
CONTINUE Dee eae dae 
MK( 35 E-¥( 95 

3005 I=1 
oxen 2) ACE) V ae I)-C4(1)*Y( J5-1) 


Bf ve) seRSLMEGS >} THEN 
c( 35) 
Bl IS) =VE( JS) +VP( IS) 


EN 
Bo. 120 I=1,N2 
)=0. 


98 


its 


ae 


3001 


LOZS 


CON 
CAnE TRANS(NG 4B M art F1,ET) 
CALL FAMIL ee 

ONG RESULT MQM Uae) 

O ged N2 


My = O06 
CONTINUE. 
GAnE ie N21, Le Use ul uy) 
C70 J) =0 
DO 3002) 1=ine 
Crea =C7(3}+U4( T, M)*U4(1I,M) 
G8 3)=8 RT(C7(J)) 
CONTINUE 
DO 8888 I=1,N1 
ei) =or 


Cc 
CONTINU 

DO 8111 I=1 

WRITE (8, — C8(I) 

CONTINUE 

FORMAT( 10X 6) 

WRITE CLE Eh SYS. OUTPUT','MOD. OUTPUT' 
FORMAT( 13X,A11,6X,A11) 

DO 976 I=1.100 

WRITE (6,8002) ¥(I),YM(I 

WRITE Vil)’ ym I 
ORMAT 6,6X,f£15.6 


,L8 
WRITE( 8,444 
WRITE( 8/4442) H1(I 
WRITE( 8,4442) H2(I 
WRITE(8,4442) H3{ I 
WRITE( 8/4442) H4/(I 
WRITE( 8, 4442) H5(I 
FORMAT le 
FORMAT( 10X/F15. 6) 
CONTI 


oe) 


O~Id elelele) QNgnaAa Qa 
(ete) 


QNQAA 


wun 
WNrH 


QnQaaAgQ 


* THIS PART PLOTS THE NECESSARY DATA 

* THAT PROGRAM CALCULATES 

KRKKKKRKKEEK pee ee es 
CAbeevPLOT (KD, SE,ba,2, F1NE 
CALL VIPLOT(KD,C8,L8,1,'TIME (, 


CALL V2PLOT( KD, MU,L8,1,'TIME : 
CALL V3PLOT( KD.V,L8,1,'TIME ' 
CALL V4PLOT( KD, 5d, nit 


Oe ti Tn tT 
-_ > 
a a a 


1 
aL 
CALL V4PLOT( KD,H4,L8,1, 
CALL V4PLOT( KD,H5,L8,1,' 


’ 
‘ 
‘ 
4 
/ 


EN 
KEEKEKRKRERRERKRKKKRRKRKEKERKRRRKRKERRKRKRRKRKRKKRRKKKE 


* THIS SUBROUTINE PRINTS THE ARRAYS AS AN* 
x MATRIX FORM * 
See Ses Fe ee Ra RK 

SUBROUTINE RESULT(H,0,P,S,T) 

INTEGER H, 0, P, S 

REARS EC 0 | 6 

PO WRITE( 8 70 T(z See) 

PRINT x(%6 Bese 25 0} 
FOL 3K. F12.6)) 

ONT INU 
RET TURN 


END 
KREEKEKKEKRERRRRKKRREKRKRRRRRRKRRKRKRKRRRKRKRRKRRKKRKE 


* THIS SUB. CALCULATES THE MULTIPLICATION* 
* TWO MATRICES * 
Soe CORE OM ORCC MC ROARED CORRE EEE Re RE EE Re ee 
SUBROUTINE CALCUL U Ve, 2k eT WO) 
INTEGER U,V 
REAL 22(U, Vie Tae, 45" Q( Y,V) 
DO 93 U 


Gi 
“Dor gt 2ko1,¥ Z,X)+W(Z,ZX)*O( ZX,X 
aie eae ) ( )+W( eel ) 


CONTINUE 
CONTINUE 
RETURN 


END 
KREERKEKRRRRKRRRRKKRKKRRRRRRERRERRKEKRKRKRERRRRRREK 


* THIS SUBROUTINE WRITES THE MATRIX NAME if 


* AND ROW,COLUMN NUMBERS. 
tek RK KKK RRR RRR KR RRR RAK RRR RR RRR RE 


100 


2-2 2 eee = ew 


i 


ond won WON 
mn 


QANAAN 


43 


Q2QaAQ 


SUBROUTINE TCOLUHN COLUMN , NAME ) 
eee Ne eer 


CHARAC NAM 
WRITE(B, 4) NATRIX' NAME 

PRINT 94 WAFRIX NAME 

FORMAT(T1,'0' 10%, gud? ,AS,T19,1X,T20,A2) 
WRITE( 8, A = ROW NUM BER. O 

FORMAT m4 Of name LOX, 713 Al}, 123, 8X,T31,12) 
WRITE( 8 NUMBER: OLUMN 

PRINT Of seu NUMBE CotuMN 

Pele, 2 Ore Ti, aia: T26,5X,7T31,12) 


END” 
KRKEEKKEKRKEKERKKKRKEKEERREKEKKKKKEKRERERREKRKKKKKKRKEKEEK 


* THIS SUBROUTINE CALCULATES THE SUMMATION * 
* OF THE SAME COLUMN ELEMENTS = 
KKK KERR RK KR KERR KR ERR RRR RRR RR KK RRR KKK KEKKKKKKKKK 


SUE ees Fee ita I2,H1,H2, Giz) 


INTEGER 
REAL ee 578i, 12) 
#30 42 H1= 
w2(1, HD )=W2(1, H2 )+ABS(G1(H1,H2) ) 
CONTINUE 
CONTINUE 
RETURN 


END 
KKK KR EKER RK RK KK RRR RK RK RE KEK KEKE RRR KEKE KEKE 


* THIS SUBROUTINE MULTIPLIES THE MATRIX * 
x BY SCALAR NUMBER * 
Frege Yc GTI TG eg TDs hs pT eT GT a la Ca gy cy yt otal GEE ct ay ah oly ch ay tha oly 
SUBROUTINE SCALAR( K1, K2,G11,,62,P1, Ga9e4) 
INTEGER Kl, G2 
REAL 3 82: ka), ie K2),P1 
DO 43 


oT en 
Ga Gil, 2) *43(c11, G2)*P1 
CONT INUE 
CONTINUE 
RETURN 


END 
KKK KEE KEKE KEK RK KKK KEKE KR RK RRR EKER RRR KK RK RK 


* THIS SUBROUTINE CALCULATES THE SUM OF * 
* THE TWO MATRICES. * 
KRRERRRRKRRRRRRRRKRRKEKRKRKKERRRKERRRKKKEKRRRREREERE 
SUBROUTINE SUM( K3,K4,15,K5,X1,X2,X3) 
INTEGER K3,K4 
REAL KICKS , ahah eke. K4) ,X3(K3,K4) 


101 


Pb 
(07 


QQaQ 


iD 
6) 


Qa|ngAaANn 


ies 
UO) 


QaagAaAaA 


46 
45 


Cc 
Cc 
Cc 


DO 46 I5=1 


X3(15,K5)= =KTC IS, K5)+X2(15,K5) 


CONTINUE 
CONTINUE 
RETURN 


END 
KKK KKK KKK KKK RRR RR KEK KERR RK RAKE KK RR RRR EEK 


* THIS SUBROUTINE COPPIES THE MATRICES. * 
KKK RK KKK KRKEE 
SUBROUTINE COPYM( K3,K4,15,K5,X1,X2) 

INTEGER K3,K4 K 
REAL X1(K3 5 Arc K4) 


DO 46’ 15=1,K3 
X1( 15, Ke) aKOCIS, KS) 
CONTINUE 
CONTINUE 
RETURN 


END 
KRHRAREKEKEKRKRKRRRERKEKEKRRRRRARKRRRRKRRKRKRKRKRKRRKRRKRKRKRKRKRREE 


* THIS SUBROUTINE FINDS THE TRANSPOSE. * 
* OF THE MATRIX. * 
Ce de de ee ee eee ee ee ee ee ee ee ed 
SUBROUTINE TRANS( K3 , K4,15,K5,X1,X2) 
INTEGER K3,K4 K 
REAL K1(K3 5 Soa K3) 


46'T5=1 
XO(KS, 12 )skicIs, K5) 
CONTINUE 
CONTINUE 
RETURN 


Pepe ges rae ae ac SPOR ea RE re Oe 
* THIS SUBROUTINE CALCULATES THE SUBTRACTION. 


* OF THE MATRICES. 


ress Sap ed ae OL Sg 


SUBROUTINE Dee Kael , KS, el, X27 43) 


INTEGER K3, 
REAL K1(K3 5 jet, K4) ,X3(K3,K4) 

DO 46’ I5=1,K3 

X3( 15, KB) ate rs, K5)-X2(15,K5) 
CONTINUE 
CONTINUE 
RETURN 


RREKKKKRKKKKRKRKRKRKRKRKRKRRKRKRKRKRKRKRKKK RRR KKKEEE 


* THs sve euets) THE ARRAYS USING - 
* THE DISSPLA PROGRAM. 2 


LOZ 


ec KREEKKKKRKEKEEKKKRKERKEKRKRKRKRKKKKKKRKRKRKRKRKKRKRKKKKKKEK 


SUBROUTINE VPLOT( T,U,M,N, LABELX, LABELY) 
CEKKKKKKEKER 


BEAL UL 800; 4) g3( 200) { 300) 7H 4), 2N( 4) 


c 
CALL AMAX(T,M, TMAX, PMAX 
CALL AMIN( TM; TMIN, PMIN 
DO 20 J = 
Yow IC) _ it 
Xi =U( 1) 
NKe CONTINUE 
Ses M, YMAX, PMAX) 
CALL AMIN( Y, M,YMIN, PMIN) 
Rite =YMIN 
Be WRIT E( 6, - YMAX , YMIN 
Weeeet on alee} 
WRITE 6. ale) 
CALL PENS fase diay 
Ge toe 6 q oMtne 
CALL (6, 1) aA 
CAL Pacer iA. 2), 
CALL 
Cc CALL NOB 
CALL AREAZD( 9 
CALL XNAME ABELK? 6 
CALL YNAME eee 6 
CALL HEADIN ( PLANT AND MODEL OUTPUTS $',100,1.,1) 
CAEL GRAE( TMIN, 'SCALE' , TMAX, YMIN, 'SCALE' , YMAX) 
GALL Spur 
DO 40 J =a ae 
a wert =u 7 0) 
30 cons thor 
ais EO. 2 oo CHNDOT 
RV M,0O) 
40 
CALL ENBPL 0) 
c Cunt Bene 
11 FORMAT Flor Acai 
13 FORMAT( " ' WOx 4612 57 
RETU 
c 


SUBROUTINE V1PLOT(T,U,M,N, LABELX, LABELY) 


ILO)! 


CERKKKKKKKKKKE 


ENEEGER PMAX)PMIN MN, LABELK, LABELY 


c 
CALL AMAX(T,M,TMAX, ey 
CALL AMIN T M; TMIN, PMIN 
DO 20 ae 
a =u(T ae 
10 Sieom 
CALL AMAX( Y,M, YMAX, PMAX) 
¥X( 5) = 
CALL AMIN(Y, M,YMIN, PMIN) 
Rie =YMIN 
a WRIT E( 6, a YMAX , YMIN 
WRITE, Be 13 TE =I. 4 
WRITE 6. 13 
CALL PANY aia 
WRITE, B. 7 IN eM 
CALL TEK 9 AX 
CALL BLOWUP 1. " 
ChErT PAGCH( 11.6. ¥ 
C CALL NOB 
CALL NaueT LR 
CALL XNAME ABEL 6 
CALL YNAME( LABELY 
CALL HEADIN ( ELE AM TER ERROR $',100,1. ,1) 
CALL ciNeAR 'SCALE' , TMAX, YMIN, 'SCALE' , YMAX) 
CALL 
DO 40 J =1, 
DO 30 I =u . 
30 CONTINUE 
J. BO. 2) ee CALL CHNDOT 
RV 0) 
40 CONTTNG 
CALL RNBPL 0) 
c CALL Donley 
11 FORMAT POMC I2. 5/7 
13. FORMAT(' ')10X/4G12. 5/ 
RETURN 
A END 
SUBROUTINE V2PLOT(T,U,M,N, LABELX, LABELY ) 
CEKKKKKEKEKE 


SAAB CEA°ORAD) cM GW) «HL S00) CHAM IM 


104 


CALL AMAX(T,M, TMAX, PMAX 
CALL AMIN( T/M, TMIN, PMIN 
0 20 J =1,N 
© 10 I’=1,M 
XI SIU ia) 
10 CONTINUE 
CALL AMAX( Y,M,YMAX, PMAX) 
CAnE AMIN( Y, M, YMIN, PMIN) 
- onfiBetE( 6rd 2 YMAX, YMIN 
WREtEY Be pis 
WRITE G4 13 
GREE A AU ‘halk 
WRITE } NMA YM IN 
Piet: bert 618 
CALL PaCet Li. 2 2). 
CALL PAGE( 1 
Cc CALL 
CaLr Ms 
CALE jReae ABED SC 6 
CALL YNAME( LABELY 
CALL HEADIN ( aE PREDICTION ERROR $',100,1.,1) 
CALL GRAF(TMIN, 'SCALE' , TMAX, YMIN, 'SCALE' , YMAX) 
CALL See 
Dortorn ave ae 
DO 30 I’=1 
Y(I =u, a) 
30 CONT 
TF( JeEO. 2) CALU CHNDOrT 
GALL CURVE (1/2 Mon 
40 oO E 
CALL ENDPL( 0) 
c CALL Roney 
11 FORMAT lope ea 2.5/ 
13 FORMAE(' ", 10%, 461215, 
RETURN 
i END 
SUBROUTINE V3PLOT( T,U,M,N, LABELX, LABELY) 
CEKKRKKEKKEKEE 
REAL U 30 0,4 800), Y¥( 800 4),YN(4) 
“ REAL U(800 4 ) bMthOM d at RBenk. x BEL 


CALL AMIN TM TMAX , ‘BMINY 
CALL AMIN( T,M,TMIN, PMIN 


IONS) 


10 CONT 
CALL AMAX( Y,M, YMAX, PMAX) 
YX(J) = YMA 
CALL AMIN( ¥,M, YMIN, PMIN) 


es B( 6,11) YMAX, YMIN 


20 
ae ch Lie 1/4} 
CALE MAX bx, bu 
CALL PMIN 
WRI oMtN 
(6/12 
CALL Ee, 2), 
Onur ROSEDA 
Cc CALL NOBR 
CALL AREA rs 
IMIG tie ABEL, 6 
CALL YNAME( LABELY 
CALL HEADIN ( Eres Aly INPUt.S 100 oie 1 ) 
CALL sshike 'SCALE', TMAX, YMIN, 'SCALE' , YMAX) 
CALL SELI 
DO 40 Le 
DO 30 I‘= 
ae: ae a 
30 CONT N 
ai EO. 23 Ge CALL CHNDOT 
RV 0) 
40 CON 
BSNL "ENDEL( 0) 
Cc CALL Roney 
11 FORMAT liek 261257 
13. FORMAT (10x, 4612.57 
RETURN 


END 
SUBROUTINE V4PLOT( T,U,M,N, LABELX, LABELY ) 
CEKKKKEKEKEE 


BRE UCORa SAL8OR) LBQB) TRL SLUM 


CALL AMAX( T,M, TMAX, PMAX 
BO 20°g St, fo fh), EM 


DO 10 ELM 
10 “ore a! 


Cc 


106 


aqaaqagqagagaAAaANd Aa 


Q 


20 


30 


40 


eo 
WH 


ial ee M, YMAX, PMAX) 
CALL pie YOM, LIN erin) 


Fe e YMAX, YMIN 
CONTINUE 


WRITE Ce? Ts ae 

ee 6. ar ,4 
SYA 
PMIN 


CALE ae YN, N. ZEN 


WRITE(6,11) ¥ MIN 
CALL TE 

CaEn Pacer is. 8 2). 

CALL 

CALL 


Gann en 


CALL YNAME tes. 
CALL HEA ONEUT TO THE PLANT Ss () 100m ies) 
CALL GROSS 
CALL GRAF(TMIN, 'SCALE',TMAX, YMIN, 'SCALE' , YMAX) 

LL SPLINE 
BO. 405 e100 

isfomc Oma = 

ith =U 
CONT NUE 
EO. 2) co FENoo 

CONT T 


SUBROUTINE AMAX( Y,N, YMAX, PMAX) 


KRREEKKKKRKRERKRKRKKKKEKKKRKRKRKRRKRRKRKRKRRKRKKRKRKRKRKRKAKKEK 


* SUBROUTINE TO COMPUTE THE LARGE ELELEMENT 


* IN A GIVEN ROW OR COLUMN IN ARRAY 
ef eS ee ee ee ee ee ee eee ee eS ee 


**x* VARIABLE DECLARATION *+= 
REAL AMAX 

INTEGER PM 

DIMENSION PCN) 


YMAX=Y oy 
DO 510 =2,N 


Og, 


ease eee LE. YMAX) GO TO 5000 
gC 


5000 CONTINUE 
5100 CONTINUE 
c MAXROW=POS 
RETURN 
END 
G 
SUBROUTINE AMIN( Y,N, YMIN, PMIN) 
S KREEKKHKHKRKKHHKEKEKKHKHHKKKEKHKKEKKHKEKKEKKKEKKKEKEKKKKEKEER 
GC * SUBROUTINE TO COMPUTE THE SMALL ELEMENT 
Cc * IN A GIVEN ROW OR COLUMN IN ARRAY 
Cc didh dy Tha th dy dhdh dh dh Th ahh Cae Uh, Te Ms th AE AE Te AD th th th Ea oh cto ty tty Pgs gt eg RO 
CG 
CG *** VARIABLE DECLARATION *** 
CG REAL AMI 
RERECER DMIN 
: DIMENSION ¥(N) 
YMINEY(1) _ 
Do 6100 1=2,.N 
TE(¥(1).GE, YMIN) GO TO 6000 
N=Y( I) 
PMIN=I 
6000 CONTINUE 
6100 CONTINUE 
ETURN 
END 


108 


APPENDIA TD 


COMPUTER PROGRAM RCIOP 


& Ke KK KKK KKRKKKKE 

CG * THESIS PROGRAM * 

GC * INDIRECT ADAPTIVE CONTROL * 

C * PARAMETER ESTIMATION * 

C * USING PROJECTION ALGORITHM * 

ce Kem KKK KKKKRKRKKKE 

C DEFINITION OF VARIABLES: 
INTEGER N1,N,L, K, ib2 J,N3,N2,31,J2,33,Mieme 
INTEGER N4‘N5,N baton 5,5; Jz. J8, 53, LS, MG, L6,a 


INTEGER DEG1, beese DEG3 sk3: 'M7, 


REAL vat i} £(p26) 1a i ath i, Pats 4) 


REAL AB : 
REAL TSP,TS TSN ie Va an 
REAL A(26,20),H 126 =a B0) 0), 36 
REAL Di 10) /T( 20 29 6 f58 6 C3 10 
REAL U3(4,1 CS 46 oie a 0 aL) Ale ) 
REAL VC( 800), 800) , ihes 
REAL SUM2,T ip. eM ) 22 50) 
REAL CON1,CON2, /SIGMALA an ee { : 266 }5 
PI=3. 141592654' 
DO 20 I =1,200 

DO 8 J soe 

SE(1,J)=0. 0 
8 continue 


M= 
WRITE( 6,701) 

701 FORMAT(‘'1'| ENTER THE ORDER OF THE SYSTEM’) 
READ(S,*) N 


L8= 2 
WRITE( 6 i292 
2301 FORMA zC NTER CONSTANT C (PROJECTION ALGORITM)') 


RET ECS. 3309 
2302 ae ov ose CONSTANT A (PROJECTION ALGORITHM) ' ) 
WRIT 


2002 FORMA c hue +L TER PLANT NUMERATOR POLYNOMIAL 


09 


2003 


2004 


Zool 
2005 


2006 
2007 


2009 
2008 


LZ 


243 


124 


2010 


Zaiod 1 


CORFRICIENTS ) 


Pa 0 

FORHAT( ey ? ty ASCENDING ORDER OF Z') 
Do 2001 I I=i,N 
eT 6, 2904) 57 
FORMAT( ‘ oan CF 2ts@ 2 =") 
continue Pe ymCS( I 

WRITE O 

RMere ! 49 Seen aa DENOMINATOR POLYNOMIAL 


COEFF INTER TS 
WRITE( 6, xt 
N ASCENDING ORDER OF Z') 


WRITE 
gaa RIGHEST COEFF. SHOULD BE l. ') 
bo “290 
WRITE 6,200 
THE 9 COREE. OF Z**(' a) 
READ( 5, ,*) C(I 
DO.112-L=1 
,vP(L)= ssn GSE1/2) 51m HELA} 8 ee 
bari ee Per o)tSint i *pi77 jtoim L*PiyZ 
DO 243 r= i 209 
ver 1, )=160. 0 
CONTINU 


DO ian 1 
wi) =O} 
U( L)=VP(L) 
CONTINUE 
Y(N1 = =ON0 
10 I= 
Y(N1 IY NI) - -C6(1)*Y(1)+C5(1)*U( I) 


CONTINU 
DO 2011 I=1,N 


O 127 t=1,N2 


dua] 


LZ3 


702 
703 
704 


705 
157 
706 


207 
708 


709 
710 


DO 122 K= N2 
TE(T, a = THEN 
ies 
P2(1,K)=0.0 
Eb te) 
CONTINUE 
CONTINUE 
DOm1121 =a No 
DO 1132 k=1 


eC 1 “&) * CHEN 
nee! K)=CON2 


PAT K ee =000 


CONTINUE 

CONTINUE 

CALL FAMILY(N2,N2,'P2' Me 
CALL RESULT(N2/N2‘1,K Za 
CALL COPYM(N2,N2,1,1L, AB, 
DOr 123. Ta 


WRITE( 6 
FORHAT( ' ? eeneey a IN DESCENDING ORDER OF Z2') 
AY ! ie HIGHEST COEFF. MUST BE 1.') 
FORHAL 79 
(= ‘*REMARK: ae SHOULD BE STABLE POLYNOMIAL’ ) 
seen 6,705) N- . . 
FOR COEFF. OF 2e*(' eye) 
ati : €3(1) 
: BAIA 7 
CG | ENTER PSTAR(Z) IN ASCENDING ORDER 
WRITE 
a | BIGHEST COEFF. MUST BE 1.') 
WRITES ee 
*REMARK-1: PSTAR(Z) IS THE DENOMINATOR 
POLYNOMIAL’ ) 
WRITE( 6,709 
FORMA OF MODEL' ) 
WRITE(S , 7) 
ce She Ne -2:PSTAR(Z) SHOULD BE STABLE 
copie 
DO 158. 121) 


WRITE( 6. 711)N-1 


Te 


Wade 
158 


241 


a7 


148 


ea 


Za? 


vonlal 


4001 


7007 


DO 


FORMAT(' ' {eee OF Z**('! 
PHU’ Pub) oe eat 


DO 241 I= 
2B TS 


114 J= 
CARE CNLGuE NOMME NS lo pale Valeo | 
CALEMCALCUR( M M)IN2. 1,L,K, V2) em, V1) ) 
TSP=CON1+V2 eRe) 

iM N2m, LL, K,V4,6r, Ul) 


CALPE ICARCUR( NOMMEN2 1,L,K 
CALL SCALAR( N2.M.1,X, 
CALL SUM(N2,M,1,K,U1,V5,0 
DO 117 I=1, N2 | 


et MO: 
NTINU 
ENG COPYM( N2 , M, I, a U} U2) 
CALL FAMIL p2 
CALL RESULT Nz, NS! I,K,P2) 


DO 7011 I=1,N 
M 


BARE} =3: 8 
BO" Oe gst )20¢ -C4( I) 

INU 
Bo 72007 | I= L, _N 

ats 
92(1)= 93(M8) 

CONT INU 
DEG1=N+1 
DEG2=2*N 
DEG3=DEG2+ 
DO 6001 I= 


2 
1,DEG3 
R(1)=0. 0 


1a 


ne) 


9003 


PU 
are 


on) 
He 


70 


714 


80 


DO 6002 K=1,DEG2 
ae I-K). LT. 0. 0). ORM Grn bac inn 


GONE RL sR Seah (Ker) #O2( Tec Re1)) 


Bo aoe TLL Ne 


ax Oo -) THEN 


END 
Ea &. GT. 0. ). AND. (K.LT.M1)) THEN 
end! 24) =D K) 


EF 

CONTINUE 

CONTINUE 

DO 61 J4=1,N 
DOr yi 
ae 

K= 

IF 


=1,N4 
=n 
1-34-11 
(es 
END IF 
CONTINUE 
CONTINUE 
DO 10 COLUM=1,N4 
SUM2=0. 


DO 70 J4=COLUM 


SUM2= guma+( AC 34, COLUM) ) **2 
CONT INU 


STGMA™SGR T( SUM2) 
ASIGMA=ABS( SIGMA MA) 
IF (ASIGMA. ut, 0001) THEN 


WRITELe, 74 
MATRIX A IS SINGULAR’) 


1 
+ 

0.) AN eee ie sis 1) THEN 
I,N+1+J4)=B(K 


END IF 

DO 80 eee ae 
G(J4)=0. 

CONTINUE 

G(COLUM)=SIGMA 


ves! 


SUM2=0. 
DO 90 J4=COLUM, N4 
SUM2=SUM2+(A(J4, COLUM) -G( J4))**2 
90 CONTINUE 
DO 100 J4=COLUM,N4 
DO 110 K=COLUM, N4 
IF (ABS(SUM2). LT. 0. 00001) THEN 
EMP=0. 


E 
ano =¢( J4 ne Core) 
ion 2B ee aan G( J4))*( A( K, COLUM) 


E 
1g Ee ee THEN 
EMP=1. +TEMP 


H( J4,K)=TEMP 
CONTINUE 
CONTINUE Z 
DO 131 J4=COLUM,N4 
DO 130 K=COLUM, N4 
TEMP=0 


DO 140 L=COLUM,N4 
TEMP=TEMP+H( J4,L)*A(L,K) 
140 CONTINUE 
C( J4,K)=TEMP 
130 CONTINUE 
131 CONTINUE 
DO 150 J4=COLUM,N4 
TEMP=0. 
DO 160 L=COLUM,N4_. 
TEMP=TEMP+H( J4,L)*E(L) 
160 CONTINUE 
T( J4)=TEMP 
150 CONTINUE 
DO 170 J4=COLUM, N4 
DO 180 K=COLUM,N4 
MG, Ky=C( 4,8) 


me 
OO 


180 C aby 
F(J4)=T(J4) 
70 CONTINUE 
70 SONA} SeCnay sac na ,v4) 
I=N4-1 
eu2 SUM2=0 
M2=I+1 
DO 30 J4= 


M2 ,N4 
SUM2=SUM2+A(1,J4)*X(J4) 
XC SCECL SUM2)/A(I,I 
Tf abst x(T)) Ed” Of 00001) THEN 


114 


141 


142 


31010 


3003 


ie 


2508 


3004 


3005 


120 


X( 1)=0. 

END IF 

f=le1 

TE( 1 chee GO TO 812 

m4 a3t1h= xT) 

CONTI be 

DO 142 I=N3,N4 
N5=I-N1 
S1(N5)=X(I) 

Oo ais 


DO! 30% Tai N 
UPC J )=UP J)+(S2(L7)-C3( 1) )*U( J-1)+S1(1)*x( gon 
PEST TL) SVC I=L) 


CONTINUE 
ee UP(J)4V( 5) ) 7 1-s2ry) 
Y fe) 

Bo” aay N 

conn yee Yl JS) -C6( 1) *¥( JS-M6) +C5( I) #U( J5-M6) 
DO 189 I=1,N1 

¥M( 1)=0.0 
CONTINU 


3 2285 —On) 
DO 2508 I=1,N 
+ 


ee TMC IS) =YM IS) - -C4(1)*YM( J5-1)+B1(1I)*V(J5-I 
CONTINUE ( Ea 
MK( 35 =-¥(95 

3005 I=1 
» HS( IS J=NKC IS )#B( 1) *V( J5- I)-C4(1)*¥( J5<1) 


mu, 35 ESTOS, 
Helga) RSL 6) THEN 


i 
5 5 
Se) = =VC( J5)+VP( JS) 


dels 


os 


aes) 


3001 


ZS 


J2= aa I 


atte jm Ot 32) 


CONTI 
DO Te” Te N2 

ET(M Td 
CONTINU 
CALL PA eCy aot M : “Ae gal ake 
CALL FAMIL ay 
eylovny RIMSIOROE NE 1) 
DO 3001 I=1,N2 

U4( I 11 N2" 


NU 
ee ‘a NZ in, U3,U1, Ua) 
j3802° 421 

© E7(Eg)=C7( 16) +041, M)*U4(1,M) 
conf SORT(C7(L6)) 
CONTINU 
WRITE(S 28% s PARAMETER ERROR' 


FORMA 
DO 8882 I= 
WRITE(8 $383 cg I) 
EORMA (15x, F12. 
CONTIN 
WRITE ok A "SYS. OUTPUT','MOD. OUTPUT’ 
FORMA 14x, Wad) 
Pere. 5. 110 Y(1I),YM(I 
WRITE (8,8002) Y(I)/YMI 
FORMAT( 9X, F10. 6,6X, £10. 6 
CONTINUE 
De 1025 I=1,L8 
SE(I,1)=¥ r) 
SE( 1, 2)=YmM( 1) 
KD I= 
CONTINU 
C8(36)=0. 65827 
C8( 37)=0. 65827 
C8( 38)=0. 65827 
C8 39)=0. 65827 
cst 40)=0. 65827 
CALL VPLOT(KD,SE,L8,2,'TIME_ ',' ») 
CALL VIPLOD( KD, CS,/L8 i, ‘TIME - 
CALL V2PLOT KD, MOS TIME 7 | 
GAteevePLorT Vlei eeu a 
CALL V4PLOT RD’ uy! ey STOR : 





LIES 


91 
22 
93 


LOE 


END 
*THIS SUBROUTINE PRINTS THE ARRAYS AS AN MATRIX FORM. * 
SUBROUTINE RESULT(H,0,P,S,T) 
INTEGER H , 0, P,S 
REA Ta aE © 
EG, 470 T(B,S),S=1 
PRINT Hees S29 


CONTINUE. 
aE 


ND 
kee THTS SUBROUTINE CALCULATES THE MULTIPLICATION OF 
TWO MATRICES. 
SUBROUTINE “caLCUL( Y V,N~Z) x Zk Zee) 
REAL 20(0 Wy, aCe *o% Q(Y,V) 
DO 93 


DO 93 sel 
20( 2 2) = _9. 9 
DOLQIeZ 
conn eye) = Zor Z, X)+W( Z,ZX) *Q( ZX,X) 
CONTINUE 


CONTINUE 
RETURN 


95°) 
FOC 4k F12.6)) 


END 
Ce **THis sUsROULINE WRITES THE MATRIX NAME AND ROW, 


wd oN won 
Ov 104) 


QQ 


: 


COLUMN NUMBERS. 
SUBROUTINE FAMIL Oe a! COLUMN , NAME ) 
INTEGER ROW, NAME 


CHARACTER*2 
WRITE(8, Har os (NAME 

PRINT 94 RI 

FORMAT( Ti TS, a Nat B12 Ee als. 1X, T200A2) 
WRTTE S82 "864° 48 NUMBER, - 

FORMAT( Ti, at a tae Pal Bay, 123 Busi 12) 
WRITE(8,96) ‘COLUMN NUMBER. ') COLUMN’ 

PRINT 96, "COLUMN NUMBER, OMN 

EFORMAT( TI, O° 72 Tome ro hos T26 Sh) sire) 


END 
***THIS SUBROUTINE CALCULATES THE SUMMATION OF THE 
absolu alues OF THE SAME COLUMN ELEMENTS. *** 
SUBROUTINE ALT 11 128 HT H2,G1,W2) 
INTEGER 11,12, 2 
REAL W2(1, £25/61(i1, 12) 


‘ley 


44 
43 


46 
45 


46 
45 


46 
45 


pio) ay) shisha 
W2(1,H2)=W2(1,H2)+ABS(G1(H1,HZ2) ) 


CONTINUE 
CONTINUE 
RETURN 


END 
***THIS SUBROUTINE MULTIPLIES THE MATRIX BY SCALAR 
SUBROUTINE SCALAR( K1,K2,G11,G2,P1,G3,G4) 


TNDEGEROKIOK2 GC 
REAL G3(K1 K2), 


DO 43 G2=1, 
DO 44'C11=1,K1 
G4(Giliy Ga=cs¢Gil ,C2)*P1 


CONTINUE 
CONTINUE 
RETURN 


Tie 
G42) Pl 


END 2 
***xTHIS SUBROUTINE CALCULATES THE SUM OF THE TWO 


MATRICES. 
SUBROUTINE 
INTEGER K3, 


K4 


SUM( K3 ,K4, 15, K5,X1,X2,X3) 


ERA Silene Ka eke K4) X3(K3,K4 
is 45 peat a ) ( ) 


CONTINUE 
CONTINUE 
RETURN 


DO 
KOris, kK 


S= ia 
5 )=K1(15,K5)+X2(15,K5) 


END 
<< THtomSUPrOU TNE CORPIES THE MATRICES. *** 
SUBROUTINE COPYM( KS, K4,15,K5,X1,X2) 


NTEGER K3,K4,15,K 
K4) MZOKs) K4) 
5 K5=1;K4 


REAL X1(K3 ) 


CONTINUE 
CONTINUE 
RETURN 


DO 46° 1I5=1,K3 
X1( 15,K5)=XK2(15,KS5) 


END 
*THIS SUBROUTINE FINDS 
SUBROUTINE TRANS(K3,K4,15,K5,X1, X2) 


INTEGER K3, 


K4,I 


THE LRANSPOoe OR ZHE MATRIX. * 


5.K 
X2(K4,K3) 


5 
REAL X1(K3,K4 
is 45 pekel Ke 


CONTINUE 
CONTINUE 


DO 
FOI 


46 I 


51 


Bae 
5)=X1(15,K5) 


JT) 


RETURN 


END 
C ***THIS SUBROUTINE CALCULATES THE SUB. OF THE 
TWO MATRICES. 
SUBROUTINE SUB, K3,K4,15,K5,xl) eee) 
INTEGER K3,K4,15 
REAL X1(K3 , A ecesit st K4) ,X3(K3,K4) 


46'15=1,K3 
R3( Ts, Ke)=K10 15, K5)-X2(15,K5) 


46 CONTINUE 
45 CONTINUE 
RETURN 


END 
kaRRRKRKRKRaRKRRRKRKRKKKRRKRKKRKRKRRRRKRKKRKRKRRRKRKKRRKKEE 


* THIS SUB. PLOTS THE ARRAYS SUsiie ie 


* THE DISSPLA PROGRAM. 
ete ee ee ee ee ee ee ee ee ee 


QaaaAaAQ 


SUBROUTINE VPLOT(1T,U,M,N, LABELA, EABEDS) 
CERKKEKREKEEE 


BEAL UL 800.4) 42( 800), ¥( 800) (4), 34) 


CALL AMIN TM, TMA, BMIN) 


Cc 


CALL AMIN( T,M, TMIN, PMIN 
DO 20 J =i, 
DO 71050, =i 


Y aut, a) 
10 CONT 
CARL AMAX( Y M, YMAX, PMAX) 
Ted j= MAK’ 
CALL AMIN Y, M, YMIN, PMIN) 


ie RITE( 6,1 11) YMAX,YMIN 


WRITE (6, ae oN a =F 2 
WRITE ge 
CALL 


CALL ri MTN EMAX 
WRITE a } YN OMAK YMIN 
CALL 8 


E 
CALL BLOWUP( 1. 2) 
CALL Eeoke 7 yon) 
BR 


20 


6 
CALL HEADIN ( PLANT AND MODEL OUTPUTS $',100,1., 


ae 


CALL Boe "SCALE',TMAX, YMIN, 'SCALE' , YMAX) 
CALL Sse link 


DOmsolrl =1,M 


¥(I =u(1,J) 
30 CONTINUE 
IF(J.£0.2) CALL CHNDOT 
CALL CORVE (T,Y,M,0) 
40 CONTINUE 
CAEL ENDEL 0) 
Cc CALL DONE j 
11 FORMAT 710%, 261255 / 
ig FoRlieme , 1OX,4G12.5/ 
RETURN 
2 END 
SUBROUTINE V1PLOT( T,U,M,N, LABELX, LABELY ) 
CRKKKKKREKE 
REAL U 30 OANA 800 Y( 800 4),YN(4 
a REAL U(80 MAX. bit N, de et epenk” Meat (s) 
CALL AMAX(T,M, TMAX, PMAX 
CALL AMIN T M,TMIN, PMIN 
DO aoa 
DO 10. =e - 
10 soci 
CALL AMAX( Y M, YMAX, PMAX) 
SO at) Nate 
ouE, AMIN Y,M,YMIN, PMIN) 
a ace 2 YMAX, YMIN 
WRITE 1,4 
WRITE alge 
CALL eae MIN chia 
WRITE “ae MAX. YMIN 
CALL 
CALL WUP(1. -2) 
ont PAGE( 1 8.5) 
Cc CALL NO 


Cann AREAZD( 9 ice 

CALL oe ABEE SC 6 

CALL YNAME( LABELY . 

CALE HEADIN ( eR EAN TER ERROR $',100,1.,1) 
CALL L GRAE( IMIN, 'SCALE' , TMAX, YMIN, 'SCALE' , YMAX) 


DO “40 J =1,N 


£20 


Dor 30 =U M 


ae } =u(I,J) 
ae CFCg EO” CALL CHNDOT 
D GORGE CT M,0) 
CONTINU 
CALL ENDPL( 0) 
“a1 FORMAT, | 10x, 2612. 5/ 
13 FORMATS ' , 10%) 2e12 5: 2%} 
RETURN 
: END 
mien SUBROUTINE V2PLOT( T,U,M,N, LABELX, LABELY) 
REAL U(800,4 800), Y( 800 4), YN( 4 
E RE REGER’P MAX, BHIN, M ds at epenk tO eety a 
CALL AMAX T, lM, TMAX, BMAX) 
CALL AMIN( TM, TMIN, PMIN 
BOO TN 
DO 10 I’=1,M 
Vi I) 20 1) 
10 CONTINUE 
Sy EEE YMAX, PMAX) 
CALL AMIN( Y, M,YMIN, PMIN) 
a ouliBetE( 6rd e YMAX , YMIN 
WRITE oe ie 
WRITE 6 613 hah 
CALL oe MEN MIN} 
WRITE( 6 1), YN OMAK Y MIN 
CAnE 
Gant ZEK Ss eD) 
CALL PAGE(11.,8.5) 
c CALL NOBRDR 
eG AREAZD( 9. 6. ) 
CATE NAMES ABET, 8} 
CALL HEADIN rade aaa PREDICTION ERROR $',100,1.,1) 
CALL GRAE( TMIN, 'SCALE' , TMAX, YMIN, 'SCALE' , YMAX) 
Colm ecrim 
Do SO SUN, M 
¥( 1) =0(1, 3) 
30 CONTINUE 


P21 


(c 


ay) 
WH 


IF(J.EQ.2) CALL CHNDOT 

CALL CURVE (T,Y,M,0) 
CONTINUE 
CALE ENDPL( 0) 
CARE Baus 
FORMAT NiO We 
FORMAT 
RETURN 
END 


SUBROUTINE V3PLOT(T,U,M,N, LABELX, LABELY ) 


AOD, He 13:27} 


CAKKKKKKKEKK 


Cc 


10 


20 


30 


40 


SBAE UC 800, 4) .2( 800),.7{ 900) (4) NV 4) 


CALL AMAX(T,M, TMAX, PMAX 
CALL AMIN T,M, TMIN, PMIN 
HO Z0 J Ne 

DO 101 


M 
ee =u(T, J) 

CONTINUE 

CAI AMAX( ¥ M, YMAX, PMAX) 

YX J YMAK ‘ 


L AMIN( YM, YMIN, PMIN) 
aa att 11) YMAX, YMIN 


ESE en oe TE L, 3} 
YHA, PA 


sbepeinyy Ie si Min 


CALL BLOWUB( 1.2 
CALL PAGE( 11. ,8. 5) 
NOB 


AREAZD( 9. 6. 
CALL oe ol ABELX, 6 
CALL YNAME( LABELY, 6 

CALL HEADIN ("EXTERNAL INPUT $',100,1.,1) 


CALL ef ee "SCALE' , TMAX, YMIN, 'SCALE' , YMAX) 
CALL SPLI 
DO 40 J zi" oN, 

DOe30 71 = 


CONTINGE =U(T, a 


ALD EO. 2 Cree woe 
a 1 CURV 0) 
CONTE 


QD 
He 
ct 
[* 
Q 


Wee 2 


1 
WH 


CALL ENDPL( 0) 


CALL DONE 

FORHAT(, ¥ nox, 2612. 27} 
FORMAT lox, 2e12) S77 
RETURN 


END 
SUBROUTINE 7) an U,M,N, LABELX, LABELY ) 


CREKKKKKKKEKE 


Cc 


10 


20 


30 


40 


=e 


WH 


REAL ue 30 800 Y( 800 4), YN( 4 
INTEGE MAR. bit N, Ys ath ae Meaty a 
CALL MENT: HY HS IS 


CALL AMIN( T,M,TMIN, PMIN 
DO.20) Tey 


CALL ie Y,M, YMAX, PMAX) 
at AMEN(Y, M, YMIN, PMIN) 
RITE(6, 1 ND Pe fistee. vaqulaeny| 


Pete acy, = i | 


Cant AMIN YN YMIN. BMAX 
WRITE( 6 yO eMbx, YMIN 


Cc E 
GALE BLOWU P 1.2) 
GALE PAGE( 11. £5 
CALL NOB 
CALL AREA2D D( 9 6.) 
CALE ae spec 8} 

HEAD iNDUT TO THE PEANT |S) 100;ie 1) 
Cane ePronm, 'SCALE', TMAX, YMIN, 'SCALE' , YMAX) 


eAur 
DO 40 J =1, N, 
DO 30 
Y y aul 1, J) 


CONTINUE 
IF(J.EO.2) CALL CHNDOT 
CALL CORVE (T,Y,M,0) 
CONTINUE 
CALL ENDPL( 0) 


CALL DONE j 
FO 710K 2612 3. 8/) 


RMAT 
FORMAT ' 10m) 4e22 57 


) 
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RETURN 
END 


c 
SUBROUTINE AMAX( Y,N, YMAX, PMAX) 
. Kee RRR KKK KK KKK KEK KEKE KRKEKEEESK 
Cc * SUBROUTINE TO COMPUTE THE LARGE ELELEMENT 
Cc * IN A GIVEN ROW OR COLUMN IN ARRAY 
e he ee ae EEN, ROW OR COLURIN IN ARRAY a bekeee 
C 
C **k VARIABLE DECLARATION *** 
c REAL AMAX 
INTEGER PMAX 
F DIMENSION Y(N) 
YMAX=Y 1) _, 
DO 5100 f= 
Be ug a YMAX) GO TO 5000 
M i) 
suis 
5000 CONTINUE 
5100 CONTINUE 
C MAXROW=POS 
RETURN 
END 
C 
SUBROUTINE AMIN( Y,N, YMIN, PMIN) 
5 keke KKK EESE 
C * SUBROUTINE TO COMPUTE THE SMALL ELEMENT 
c * IN A GIVEN ROW OR COLUMN IN A 
C Ebb hy Ty aS Phy lh SS AE ty lta ld Os En ob I hb GS IG TH ch Mth TH ca os ly eg td Oly OT 
c 
C *** VARIABLE DECLARATION ec 
C REAL AMIN 
INTEGER PMIN 
P DIMENSION Y(N) 
YMIN=Y 1) _, 
DO 610 
TFC Y( 1). {GE, YIN) GO TO 6000 
MIN=I 
6000 CONTINUE 
6100 CONTINUE 
RETURN 
END 
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10; 


dig 
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